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We study the effect of Coulomb interaction between charge carriers on the properties of graphene 
monolayer, assuming that the strength of the interaction is controlled by the dielectric permittivity 
of the substrate on which the graphene layer is placed. To this end we consider the tight-binding 
model on the hexagonal lattice coupled to the non-compact gauge field. The action of the latter 
is also discretized on the hexagonal lattice. Equilibrium ensembles of gauge field configurations 
are obtained using the Hybrid Monte-Carlo algorithm. Our numerical results indicate that at 
sufficiently strong coupling, that is, at sufficiently small substrate dielectric permittivities e < 4, 
and at sufficiently small temperatures T < 1 ■ 10 4 K the symmetry between simple sublattices of 
hexagonal lattice breaks down spontaneously and the low-frequency conductivity gradually decreases 
down to 20 — 30% of its weak-coupling value. On the other hand, in the weak-coupling regime (with 
e > 4) the conductivity practically does not depend on e and is close to the universal value oo = 1/4. 
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Introduction 

Graphene, a two-dimensional crystal with hexago- 
nal lattice formed by carbon atoms, has attracted a 
lot of attention in recent years both as a novel mate- 
rial with many unusual properties 1 ' 2 and as a unique 
laboratory which allows to study numerous quantum- 
field-theoretical phenomena in desktop-scale experiments 
(see for a review and further references). 

Electronic transport properties of graphene are of par- 
ticular interest for industrial applications. Theoretical 
considerations within the tight-binding model of crystal 
lattice show that in the absence of interactions and at low 
energies charge carriers in graphene behave as massless 
Dirac fermions with an effective "speed of light" being 
equal to the Fermi velocity vf ~ c/300 ' ' . It turns out 
that the conductivity of these massless Dirac fermions 
takes a universal value (To = 1/4 e /h in the limit of zero 
temperature and in the absence of interactions ' -12 . This 
value, however, should strongly depend on the measure- 
ment procedure and on the geometry of a sample 1 ". 

Due to the smallncss of vp, electromagnetic interac- 
tions between these Dirac fermions are well described 
by the instantaneous Coulomb potential. However, for 
the same reason the resulting effective field theory turns 
out to be strongly coupled with a coupling constant 
a = e 2 /vp ~ 300/137 ps 2. Thus one can expect that 
Coulomb interactions between electrons can significantly 
modify the properties of graphene such as the quasipar- 
ticle spectrum and the conductivity. For graphene on a 
substrate with dielectric permittivity e Coulomb inter- 
action is screened so that the coupling decreases by a 
factor ^rj. This provides a practical way to control the 
interaction strength. 

Experimental studies of the conductivity of suspended 



graphene 13,14 , for which the coupling constant is max- 
imal, suggest the existence of a gap in the quasiparti- 
cle spectrum with the width of order of lOmeV 11-14 . 
The opening of a gap due to strong Coulomb interac- 
tions for a > 1 is also supported by analytical calcu- 
lations based on the solution of the gap equation 9 ' 15,16 
and on the strong-coupling expansion in lattice gauge 
theory 1 r_19 . The transition to the gapped phase is likely 
to be of the second order 1,1 . Within the effective the- 
ory of Dirac quasiparticles the opening of a gap in the 
spectrum is accompanied by a formation of the fermion 
chiral condensate ('tpip) 6 ' 16,20 . In terms of the original 
tight-binding lattice model, such condensate corresponds 
to the difference of the charge carrier densities on the two 
simple sublattices of the hexagonal lattice 6,18,19 . 

On the other hand, more recent measurements 21 ' 22 in- 
dicate the absence of a gap in the quasiparticle spectrum 
of suspended graphene as well as logarithmic divergence 
of the Fermi velocity near the Fermi points, in agreement 
with analytical calculations based on renormalization- 
group techniques , dynamical mean- field theory" and 
expansion in the large number of fermion flavors 20 ' 25 . 

In view of such uncertainties both in the experimen- 
tal measurements and analytical calculations, it seems 
natural to turn to the first-principle numerical methods, 
such as Monte-Carlo simulations on the lattice. This line 
of research has been actively pursued recently. In 26-30 
the low-energy effective field theory of Dirac quasiparti- 
cles in graphene was studied numerically by using 2 + 1- 
dimensional staggered lattice fermions coupled to 3 + 1- 
dimensional non-compact Abelian lattice gauge field. A 
hint on the second-order semimetal-insulator phase tran- 
sition at a c = 1.11 ± 0.06 associated with the opening 
of a gap in the energy spectrum and spontaneous break- 
ing of the chiral symmetry of Dirac fermions was found. 
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In 31 ' 32 a similar model with staggered fcrmions and a 
contact interaction term instead of the Coulomb poten- 
tial was studied, and a phase transition with respect to 
the coupling constant separating the gapless conducting 
weak-coupling phase and the gapped insulating strong- 
coupling phase was also observed. A study of the finite- 
temperature phase transition in this model was reported 
in 33 . A phase transition of the Bcrezinskii-Kostcrlitz- 
Thouless type was found at T = 0.055 (2) A , where Ao 
is the width of the gap in the spectrum. 

A common strategy of the works 26-33 was to use the 
value of the fermion condensate as an order parameter. 
Recently, the results of the direct numerical measure- 
ments of the conductivity of graphene lattice effective 
field theory with staggered fermions were reported in 34 . 
It was found that the DC conductivity obtained from the 
Green-Kubo relations indeed rapidly decreases when the 
fermion condensate is formed, in agreement with theo- 
retical expectations. 

The effect of lattice artifacts of staggered fermions on 
the flavor symmetry breaking in graphene effective field 
theory was discussed in , with the conclusion that lattice 
simulations with staggered fermions might yield some- 
what lower value of the critical coupling constant than 
in the continuum graphene effective field theory Thus 
the influence of such lattice artifacts could probably shift 
the critical coupling constant of the semimetal-insulator 
phase transition above the value of the coupling constant 
in suspended graphene. Such a shift would then explain 
the fact that the insulating state of the suspended mono- 
layer graphene was not observed experimentally 21 ' 22 . 

Lattice regularization of the effective field theory of 
graphene used fri 26 ~ 35 involves two approximations: one 
first starts from the tight-binding lattice model of the 
electron transport in graphene 1 ' ' and derives the low- 
energy effective theory of Dirac fermions. This theory 
is then again discretized on the lattice by using suitable 
lattice fermions which reproduce the Dirac spectrum at 
low energies. However, since the original model is al- 
ready formulated on the hexagonal lattice, it is tempt- 
ing to circumvent these two approximations and per- 
form direct simulations of the tight-binding model with 
Coulomb interactions included. This possibility has been 
recently discussed in * . Such simulations, while tech- 
nically being even simpler than simulations with stag- 
gered fermions, have several crucial advantages. First, 
they allow to study the patterns of spontaneous symme- 
try breaking which are specific for the hexagonal lattice, 
such as the Kekule distortion. Despite the fact that sub- 
lattice symmetry is not broken in this case, a gap in the 
spectrum might develop 19 '' 9 . Second, since the lattice 
spacing is fixed in the tight-binding model, simulation 
results can be unambiguously compared with experimen- 
tal data. Finally, all the symmetries of the tight-binding 
model are explicitly preserved. As discussed in Section 2, 
the latter has a U (1) <8 U (1) symmetry associated with 
the conservation of the total numbers of charge carriers 
with different spins as well as the discrete sublattice sym- 



metry. Sublattice symmetry can be explicitly broken by 
the staggered potential to, which plays the role of the 
Dirac mass at low energies. In contrast, massive stag- 
gered fermions have only single global U (1) symmetry 
associated with total charge conservation 35 ' 40 ' 41 . Thus 
simulations of the tight-binding model are free from lat- 
tice artifacts and can serve as a completely independent 
cross-check of simulations with staggered fcrmions. In 
particular, one can estimate the influence of lattice arti- 
facts of staggered fermions on the simulation results. 

In this paper we report on the results of such direct lat- 
tice Monte-Carlo simulations of the tight-binding model 
of graphene on hexagonal lattice with Coulomb interac- 
tions. To account for the latter, we couple the tight- 
binding model to the 3 + 1-dimensional non-compact 
Abelian lattice gauge field, as in 26-35 . In contrast to 
the approach of 5 '' 38 , where the interactions are treated 
by applying the Hubbard-Stratonovich transformation, 
the resulting action is local and gauge fields outside of 
graphene plane can be efficiently sampled by a heat-bath 
algorithm 12 . As we demonstrate in Subsection 1.3, our 
discretization of electric field on the hexagonal lattice re- 
produces the continuum Coulomb potential with a very 
good precision. An additional advantage of the use of 
3 + 1-dimensional gauge fields is the possibility to study 
electromagnetic interactions of multilaycred graphene, 
for example, Casimir forces. 

We focus on the study of spontaneous breaking of sub- 
lattice symmetry as well as on the direct measurements of 
electric conductivity of graphene. In agreement with the 
results of 26-35 , we find that sublattice symmetry is spon- 
taneously broken for coupling constants a > 1 , which cor- 
responds to the substrate dielectric permittivities e < 4, 
and for sufficiently small temperatures T < 1.3 ■ 10 4 K. 
At higher temperatures sublattice symmetry is not bro- 
ken for the physical values of e, e > 1. 

It should be stressed, however, that in this paper we 
only consider the tight-binding model at finite temper- 
ature with the purpose of finding the range of lattice 
parameters which describe the low-temperature phase of 
this model. For this reason we also do not take into 
account the thermal fluctuations of the hexagonal lat- 
tice itself, which should become irrelevant for sufficiently 
low temperatures. As we will demonstrate below (see 
Subsection 1.6), realistic lattice parameters correspond 
to quite high temperatures of the electron gas of order 
of 1 eV ~ 10 4 K. Simulations at significantly lower tem- 
peratures are computationally very expensive, and it is 
important to give an upper bound for the temperatures 
which still describe the low-temperature phase in order 
to make an optimal choice of lattice parameters for prac- 
tical simulations. As we will show in Subsection 2.2, the 
critical temperature which separates the low- and the 
high-temperature phases of the tight-binding model can 
be estimated as T c w 1.3 • 10 4 K. Thus one can hope 
that for the smallest temperature which we use in our 
simulations, T = 8.8 ■ 10 3 K, the low-temperature prop- 
erties of the tight-binding model are already reproduced 
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with a sufficiently good precision. Since at temperatures 
of order of 10 K a real graphene monolayer would be 
destroyed due to thermal fluctuations of the lattice 43 , a 
detailed study of the finite-temperature phase transition 
within the tight-binding model (without taking into ac- 
count the phononic degrees of freedom) seems to be of 
purely academic interest only. 

We extract the conductivity from the correlators of 
electric current densities with the help of the Green- 
Kubo relations. In order to estimate the conductiv- 
ity in a model-independent way, we consider the low- 
frequency conductivity a smeared over frequencies w < 
T. We observe that in the strong-coupling regime the 
low-frequency conductivity quickly decreases with e down 
to 20 — 30% of its weak-coupling value at all tempera- 
tures which we have considered. On the other hand, in 
the weak-coupling regime (e > 4) the conductivity prac- 
tically docs not depend on e. 

The outline of the paper is the following. In Section 1 
we give a detailed description the geometry of hexagonal 
lattice and its extension to the (3 + l)-dimensional space, 
and describe the lattice actions for the gauge field and 
for fermions which we use in our simulations. The de- 
tails of our simulation algorithm and the choice of lattice 
parameters are also discussed in this Section. In Section 
2 we discuss spontaneous breaking of sublatticc symme- 
try. Our numerical measurements of the conductivity 
are summarized in Section 3. In Section 4 we give some 
concluding remarks on our results. Some technical de- 
tails and supplementary material (such as the calculation 
of the current-current correlator for the non-interacting 
tight-binding model) are relegated to the Appendices. 



1. LATTICE ACTION AND SIMULATION 
METHOD 

1.1. Tight-binding model of graphene with 
electromagnetic interactions 

With a good precision the electronic properties 
of graphene can be described by the tight-binding 
Hamiltonian 5 ' 7 

Htb = -K X! ( & t,x & *,y + & t,Y °><r,x) > (!) 

CT=f4 <XY> 

where summation goes over all neighboring sites X, Y 
on the hexagonal lattice with hexagon side a = 0.142 nm 
and k « 2.7 eV is the hopping energy for carbon tt or- 
bitals. a a x and a a .x are the creation and annihilation 
operators for non-relativistic electrons with spin a =^, ^. 
This Hamiltonian describes electron hopping between 
nearest neighbor atoms only. The hopping energy for 
hopping between next to nearest neighbor atoms is much 
smaller than k " and we neglect it here. 

Since in the ground state graphene is electrically neu- 
tral, there should be on average one electron per lattice 



site. We assume that the ground state of the free tight- 
binding Hamiltonian (1) is fixed by the conditions 6,3 ' 38 

a^ x |0> = 0, a\ x |0) = 0. (2) 

Thus there is one electron with spin a =4- at each lattice 
site. For the ground state fixed by (2), it is convenient to 
define the creation and annihilation operators for "parti- 
cles" and "holes" as 37 ' 38 

ipt,x = a-f, x , i>±,x = ±4,x- ( 3 ) 

In the definition of ipl,x, we take the plus sign for lattice 
sites which belong to one simple rhombic sublattice of 
the graphene hexagonal lattice and the minus sign for 
lattice sites on another simple sublattice 37 ' 38 . Now the 
ground state satisfies the standard condition ^,x |0) = 
0, i>i,x 1 0) = and the charge operator reads 

qx = V^,x it,x ~ ^\,x $i,X- (4) 

In other words, we interpret the absence of electron at 
some lattice site as the positively charged hole, and va- 
lence electrons in graphene play the role of the Dirac 
sea. Obviously, the tight-binding Hamiltonian in terms 
of the new operators has the same form as (1) with an 
additional shift of energy which docs not affect physical 
results. It should be stressed that since in the partition 
function (and hence also in lattice Monte-Carlo simula- 
tions) we anyway sum over all possible states of the the- 
ory, the particular choice of the "perturbative" vacuum 
state (2) is only a matter of convenience. 

Coulomb interaction is described by the interaction 
Hamiltonian 

Hi = \ J2e 2 /r(X,Y)q x q Y , (5) 

X,Y 

where r (X, Y) is the distance between lattice sites X 
and Y, e 2 w 1/137 is the electron charge. Through- 
out the paper we use the natural system of units with 
c = h = kg = I. A common way to simulate theo- 
ries with four-fermion interaction of the form (5) is to 
apply the Hubbard-Stratonovich transformation and to 
sample the fictitious Hubbard-Stratonovich field. For the 
long-ranged three-dimensional Coulomb potential (5) the 
resulting action is nonlocal 3 '' 38 . 

Here we adopt a different strategy and consider the 
tight-binding model (1) coupled to the real electromag- 
netic field. This coupling is introduced by using the stan- 
dard Peierls substitution within the tight-binding Hamil- 
tonian (1)36,44,45. 

a-l,x °-<r,v °t,x ex P °-<?,y, (6) 

Y 

where 6xy = —&yx = e J dx l Ai is the operator of the 

x 

integral of the electromagnetic vector potential A t along 
the lattice bond which joins the sites X and Y. 
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Furthermore, we use the operators Oxy and the mo- 
menta canonically conjugate to them to approximate the 
Hamiltonian of the electromagnetic field in continuous 
space 

Hem = Jd 3 r (e 2 (r) + (rot A (r)) ^ , (7) 

where Ei, i — 1,2,3 is the operator of the electric field 
strength, and to construct the corresponding lattice ac- 
tion. It turns out that such a discretization of the electro- 
magnetic field reproduces the continuum Coulomb poten- 
tial with a very good precision (see Fig. 4 in Subsection 
1.3). The Hamiltonian (7) should be also supplemented 
with the Gauss law constraint 

VE(r)=4ireJ2 q X 6(f,f x ). (8) 

x 




FIG. 1: Dispersion relation E ykj for the tight-binding 

Hamiltonian (9) with the staggered potential of strength m at 
different ratios m/n. The points on the plot and the inset il- 
lustrate the filling of the graphene Brillouin zone with discrete 
lattice momenta for f8 x f8 lattice with periodic boundary 
conditions. 

Let us also note that the spectrum of the tight-binding 
Hamiltonian (I) has two zero-energy states which corre- 
spond to the Fermi points (see Appendix A for a more 
detailed discussion). These zero modes might result in 
certain singularities of the effective action of electromag- 
netic field, namely, in the appearance of zero modes of 
the fcrmionic hopping operator. These zero modes make 
standard numerical simulation methods such as Hybrid 
Monte-Carlo inapplicable 26 ' 27 ' 31,42 . In order to prohibit 
the existence of zero modes, it is convenient to introduce 
the staggered potential which is equal to + m a a a x a CT ,x 
for lattice sites which belong to one simple sublattice of 
the hexagonal lattice and — m a a/ a x a a ,x for sites of the 
other sublattice. As we will see from what follows, for 
the purpose of numerical simulations it is convenient to 
take = — = m. The addition of such potential to 
the tight-binding Hamiltonian (f) opens a gap of width 
2 m in its spectrum and prohibits the existence of zero 
modes (sec Appendix A), however, at the cost of explicit 



breaking of sublattice symmetry. As we show in Ap- 
pendix B, this gap cannot be closed due to interactions, 
thus at nonzero m Monte-Carlo simulations are possible 
for any value of the coupling constant e 2 in (5). At low 
energies the coefficient m plays the role of the mass of 
Dirac quasiparticles in graphene*. In order to describe 
the physics of massless fcrmions, wc should extrapolate 
simulation results to m = 0. This situation is very sim- 
ilar to lattice QCD simulations, which are only possible 
at nonzero quark masses and the chiral limit is reached 

only by extrapolation. Dispersion relation E (k\ for the 

tight-binding Hamiltonian (1) with the staggered poten- 
tial of strength m is shown on Fig. f at different ratios 
to/ k. The points on the plot and the inset illustrate the 
filling of the graphene Brillouin zone with discrete lat- 
tice momenta on the two-dimensional toric lattice made 
of f 8 x f 8 hexagons (see Subsection f .2 and Appendix A 
for more details). 

Finally, taking into account all the refinements of 
the tight-binding model discussed above and using the 
fcrmionic operators introduced in (3), we arrive at the 
following Hamiltonian: 

H t b = -K ^ {^l-x CX P (i'^y) 4>a,Y+ 

ff=t,4.<xy> 

+i>l tY exp (±i§Yx) i>a,x) + 

+ Yl m ^°,xA°,Xi - S ^2 m ^i,x a ^,x a (9) 
cr=t,4- Xi o-=f,4- x 2 

Since particles and holes have opposite charges, in the 
first term in (9) one should take the plus sign before 
9xy and 6yx for a =t and the minus sign for a =\.. 
In the second term, summation over X\ and X2 denotes 
summation over the sites of two simple sublattices of the 
hexagonal lattice. 

A starting point for lattice Monte-Carlo simulations 
is the path integral representation of the partition func- 
tion and operator expectation values for the tight-binding 
model interacting with electromagnetic field: 

Z = Tr cxp (-H/t\ (fO) 
(Oi (n) ... O n (r«)) = 
= Z- 1 Tr (6x {n) ...6 n (t„) exp , (11) 

where H = H t b + H em , T is the temperature, the trace 
is taken over the joint Hilbert space of the fcrmions and 
the electromagnetic field, 0\ , . . . , O n are some operators 

and O (t) = exp f—r H^j O exp (t H^j . As usual, exactly 

zero temperature cannot be reached in lattice Monte- 
Carlo simulations, but a reasonably small value can be 
achieved by using lattices with sufficiently large size in 
Euclidean time direction. 

We construct the lattice approximation to the path 
integral representation of the partition function (10) in 
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Subsections 1.3 and 1.4 below. In Subsection 1.3 we ap- 
proximate the trace over the states of the electromagnetic 
field and in Subsection 1.4 - over the fermionic degrees 
of freedom. Before that, in Subsection 1.2 we describe 
in details the geometry of the hexagonal lattice and its 
extension to the three-dimensional space. 

1.2. Geometry of hexagonal lattice and its 
extension to the three-dimensional space 

In order to perform lattice Monte-Carlo simulations, 
we should somehow compactify the hexagonal lattice 
on which the tight-binding model is defined. Here we 
consider lattices which have the topology of the torus. 
An example of such hexagonal lattice which consists of 
L x x L y = 6 x 4 hexagons is shown on Fig. 2. Since 
hexagonal lattice can be thought of as a composition of 
two rhombic sublattices, it is convenient to classify the 
lattice sites which belong to these sublattices as either 
"even" or "odd" sites. All nearest neighbors of an even 
site are odd sites, and vice versa. On Fig. 2 even sites 
are marked with red circles, and odd sites - with green 
crosses. 




FIG. 2: Cartesian and rhombic coordinate axes and coor- 
dinate grids for hexagonal lattice covering the torus of size 
L x x L y — 6 x 4 in the two-dimensional Euclidean space. 

Lattice sites can be enumerated using the integer- 
valued coordinates £i = . . . L x — 1, £2 = . . . L y — 1 
which label the sites of one of sublattices, say, even sites. 
The corresponding coordinate axes and coordinate grid 
are shown on Fig. 2 with solid and dashed red lines. 
Numbers in parentheses near even lattice sites are their 
coordinates (£i,£2)- Coordinate system for odd sites is 
the same as for the even sites, but its origin is shifted 
along the hexagon edge which is perpendicular to the 
£2 axis. Altogether, we characterize each lattice site of 
the hexagonal lattice by two integer-valued coordinates 
£ = (£1, £2) and a label s = a, j3, where s = a stands 
for even lattice sites and s = /3 - for odd lattice sites. 
We also introduce the Cartesian coordinate system with 
coordinates x and y in the graphene plane, so that the X 
coincides with £1 axis of the rhombic coordinates. The 



axes of this coordinate system are shown on Fig. 2 with 
blue solid lines. Cartesian coordinates of the lattice site 
with rhombic coordinates (s, £) are 

x = \/3a£i + V3/2a^2 + V3/2aS St p, 

y = 3/2a&-l/2a5 s ,p, (12) 

where a is the lattice spacing, that is, the length of 
hexagon edge. Correspondingly, the area of the unit cell 
of the rhombic lattice in cartesian coordinates is equal to 
the hexagon area 3v / 3a 2 /2. 

In order to embed our lattice into Euclidean space with 
torus topology, we identify the opposite sides of rectan- 
gle of size V3L X x 3/2 L y in Cartesian coordinates, as 
shown on Fig. 2. Such identification implies the following 
identification of rhombic coordinates 4 ' 3 : 

(£1 + £*,&)-> (£1,6), 
(£i,£ 2 + i y )^(£i + V2,6)- (13) 

We assume that the links of our hexagonal lattice are 
always directed from even sites to odd sites, as illustrated 
on Fig. 2. Correspondingly, we label them by the coordi- 
nates £ = (£1, £2) of the even site from which they origi- 
nate and the direction number b = 0,1,2, such that the 
link with coordinates £ in direction b goes from the site 
with coordinates (a, £) to (f3, £ + Pb) (modulo the identi- 
fication (13)). Here we have introduced the following set 
of three vectors in rhombic coordinates: 

po = (0,0), Pl = (-l,l), p 2 = (-l,0). (14) 




FIG. 3: Direct product of two-dimensional hexagonal lat- 
tice with lattice spacing a and rectangular lattice with lattice 
spacing Az which fill the three-dimensional space. Dark blue 
faces which form a triangular prism are the plaquettes of the 
dual lattice. 
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In order to describe electromagnetic fields, which prop- 
agate in the three-dimensional space, we also introduce 
the coordinate z for the direction perpendicular to the 
graphene layer. In the following we assume that the lat- 
ter is situated at z = 0. We denote the vectors in the 
three-dimensional space by arrows: r = (x, y, z). We dis- 
cretize the z coordinate into the intervals of size Az, so 
that the three-dimensional space is covered by a direct 
product of the hexagonal lattice in the x, y plane and 
the regular rectangular lattice in the z direction, as illus- 
trated on Fig. 3. Shifts along the links of this lattice are 
described by the vectors 

eb = (Vs/2a,-l/2a,0j , ei = (0, a, 0) , 
e 2 = f-\/3/2a,-l/2a,o) e z = (0,0,Az). (15) 

It is also convenient to introduce the dual lattice 
with lattice sites which are situated above the centers 
of hexagons of the original lattice and which are shifted 
by Az/2 along the Z axis. The projection of the links 
of this dual lattice is shown on Fig. 2 with dashed lines. 
Now with each lattice link parallel to the graphene plane 
we can associate the rectangular plaquette of the dual 
lattice (with size ^/3a x Az) which is orthogonal to it. 
Correspondingly, each link which goes in the z direction is 
associated with some plaquette of the dual lattice which 
is parallel to the graphene plane and which has a form 
of equilateral triangle with side y/3a. Plaquettes of such 
dual lattice are also shown on Fig. 3. 



1.3. Lattice action for the electromagnetic field 

We discretize the Hamiltonian (7) on the three- 
dimensional lattice described above in Subsection 1.2. As 
discussed in Subsection 1.1, upon discretization the vec- 
tor potential is replaced by its integrals along the lattice 
links: 



b(£,z) = e / du e\ ■ A { [x {a, £,z)+u e b ) 



9 Z (s, £,z) = e j du el ■ A, (x (s, f , z) + u e z ) . (16) 
o 

From now on, we replace the abstract labels X, Y of lat- 
tice sites used in Subsection 1.1 by the coordinates s,£,z 
introduced in Subsection 1.2. Since for different s, £, z 
we take vector potential in different points, all variables 
(16) should be considered as independent. Note also that 
while the operators 9b (£, z) arc associated only with even 
lattice sites, the operators Z (s, £, z) are associated with 
both even and odd sites. 

The momentum operators canonically conjugate to 
Ob (£, z) and 9 Z (£, z) can be constructed as operators of 



electric field flux through the plaquettes p* of the dual 
lattice which are dual to the corresponding links: 



TTfc (£, Z) 



1 

Aire 



ddiE 1 ,6 = 0,1,2, 



p*±e b 

= 4^ / da ^\ (17) 

p* -Le z 

where da is the element of area on the dual plaque- 
ttes. Since by construction the operators 7Th, % z sat- 
isfy the canonical commutation relations with the op- 
erators 9b, 9 Z , they can be represented as differential op- 
erators 7T b (£, Z) = ~i gg^y 7T Z (S, t Z) = -i gg^ iZ) 

on the Hilbert space of functions of the variables 9b (£, z), 
9 Z (s,£,z). For notational convenience, let us also intro- 
duce the redundant set of field and momenta operators 
8b (s, £, z) and Wb (s, £, z) associated with each lattice site, 
either odd or even: 

9 b {a,i,z) = 6 b 
9 b (0,£,z) = -§ b (t-p b ,z) 

TT b (a,€,z) = 7T b (£,z) 

Ttb(P,S,z) = -*b(S-pb,z) (18) 

Let us first consider the discretization of the electric 
part of the Hamiltonian (7). Below we will see that 
since the charge carriers in graphene move with veloc- 
ities which are much smaller than the speed of light, the 

magnetic term ^rot Aj in the Hamiltonian can be ne- 
glected. To the leading order in a and Az we can write 

, . \f?> Az _ p,,-,, 
Kb (£, z) « — e b - E(r (a, £, z)) 



47re 
3^3a 2 
167re Az 



e z -E(r(s,£,z)), (19) 



where we have taken into account that the areas of the 
plaquettes p* dual to lattice links in the graphene plane 
and perpendicular to it are equal to y^3a Az and 3v ^ a , 
respectively. Using the identity 



V- - - 3 a I r 

2_^eb®eb = —— I 



e z (8) e z 



6=0 



(Az) 



(20) 



we can now express the square of the electric field in 
graphene plane as 



El{r{a,Z,z)) + E 2 y {r{a,t,z)) 
32tt 2 e 2 



9 a 2 Az 2 



(21) 



b=0 



The integral over the three-dimensional space in (7) can 
be also approximated by the sums over the vertices of 
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the hexagonal lattice: 



3>/3a 2 Ax 



s—ct,/3 z 



Finally, we arrive at the following discretization of the 
electric part of the Hamiltonian (7): 



1 

87 



d 6 r 



* + 



87re Az v ^ „ 2 



(23) 



s— a, ft 



By similar reasoning one can show that the discretiza- 
tion of the magnetic part of the Hamiltonian (7) should 
contain the square of some combination of the opera- 
tors 9 b (£, z) and 9 Z (s, £, z) multiplied by factors of order 

8tt e a Az aI1< ^ 8tt e a a ■ 

Let us now take the trace over the states of electro- 
magnetic field in the partition function (10). To this end 
we use the standard Feynman-Kac transformation and 
rewrite 



exp -H/T 



L T -1 

exp(-ffAr 

t/At=0 



(24) 



where At = (TL T ) , and insert the identity operators 
I = JVA(r) \A(r}) (A(r) | decomposed into the com- 
plete set of eigenvectors \A (r)) of the operators Ai be- 
tween these factors. Upon discretization (23), the Hilbert 
space of the discretizcd theory is equivalent to the space 



of all functions of link variables 9 b (£, z) and 9 Z (s, £, z), 
and the decomposition of identity reads: 



z) 



/= n (n / Mbit**)) f n / 

\9 b (£, z) , 9 Z (s, z)) (9 b z) , Z (a, £, z) | ,(25) 

where \9 b (£, z) , # z (s, £, z)) are eigenvectors of the oper- 
ators 6*;, (£, z) and a (s, £, z). For the sake of brevity, we 
will denote the integrals over 9 in (25) as T>6, and the 
corresponding eigenvectors as \6). 

The decomposition of identity (25), however, contains 
non-physical states which violate the constraint (8). In 
order to get rid of such states, one should also insert the 
projectors V on the physical Hilbert space between the 
exponents in (24). Integrating the constraint (8) over 
the unit cell of the dual lattice which encloses the lattice 
site with coordinates (s, £, z) and taking into account the 
definitions (17), we see that the discretized version of the 
constraint (8) becomes 



E^ b +7T* 



h=0 



-n z (s, £,z- Az) = q (s, £) 5 (z, 0) , 



(26) 



where <?(s,£) is the charge operator (4) at lattice site 
with coordinates (s,£) and we have taken into account 
that graphene layer is placed at z = 0. 

It is convenient to rewrite the projection operator as 
an integral over the Lagrange multiplier field <fi (s, £, r, z), 
which becomes the electric potential field in the path 
integral formalism. The matrix element of the r-th factor 
in (24) between the states (8 (r) | and \6 (r + At)) can 
be now written as 



(9 (t) I V exp (-H em At) \9 (r + At)) = JJ / # (s, £, t, z) 
' 0) I exp £ i0 (S, £, T, z) I £ 7T b (s, £, 2) + TTz (s, £, ^) - n z (s, £, z — Az) — q (s, 5 (z, 0) J X 

\6=o / / 

(27re 2 Ar . 87re 2 AzAr A ,, „ . ) ... ... 

-^ fZ7 E^(^)-^ f ^E^(^^) |«(r + Ar)> x 
£,z,6 s,£,z / 



cx p(-E-^( rotA ^^ T+Ar ' z )]) 2 )> 



s,t,z 



(27) 



where rot A [9 (t + At)] is the lattice discretization of the rotor of the vector potential and we have temporarily 
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omitted the fermionic part of the Hamiltonian H. Follow- 
ing the standard procedure, we have also approximated 
the exponent of the Hamiltonian by a product of expo- 
nents of the kinetic and potential terms. Evaluating the 



remaining matrix element which contains the exponen- 
tials of momentum operators, we arrive at the following 
expression: 



exp 



3\/3a 2 



(0(T)|Pexp(-# em Ar) |0(t + At))= [] / dct>{s,^r,z) 

% E fa ( a ' r ' z ) - ^ * + ^ r ' z ) + 9 » & r ' z ) - 0& & r + At ' z )) 2 ~ 



87re 2 A 



6,5,z 



327re 2 AzAt 



2 (0 (*, £, r, z) - </> (s, £, r, z + Az) + ^ (s, e, r, z) - 0, (s, f , r + At, z)) 2 
x exp f-i VJ , (., f , T , , - 0) , (., - £ A T '* « + ^ »)»• ) , 



r 



(28) 



Now we can collect all such factors into the path in- 
tegral over the lattice fields 6b (£, r, z), 6^ (s, £, r, z) and 
( s ) T i z )- Before writing down the final result for the 
lattice action, let us consider the proportion between the 
electric (first exponential) and the magnetic (last term in 
the second exponential) terms in the path integral weight 
(28). The coefficients before finite differences of lattice 
fields in the electric and the magnetic terms in (28) are of 
order -; — % a . and . , respectively. In lattice sim- 

ulations, we should choose At such that kAt <§C 1, but 
kAtL t ^> 1. For realistic simulations, L T ~ 10 1 , thus 
k At ~ 10 _1 . Taking into account that na = 1.946- 10 -3 
and e 2 « 1/137, we conclude that the coefficients before 
the electric and the magnetic terms in the lattice action 
are of order of 10 _1 and 10 3 , respectively. Therefore, in 
Monte-Carlo simulations the fluctuations of the spatial 
component of the gauge field 6^ (£, r, z) with nontrivial 
field strength are suppressed by two orders of magnitude 
in comparison with fluctuations of the electric potential 
4> (s, £, z, t), and one can disregard them, assuming that 
rot A [8 (£,t + At, z)] is effectively equal to zero. By a 
gauge transformation one can then set all the spatial link 
variables 6^ (£, r, z) to zero. 



Physically the dominance of the electric part of the ac- 
tion means that we adjust At so that the ratio a/ At is 
comparable with characteristic velocity of charge carriers 
in graphene vf = 3/2«;a ps 1/300. Since this velocity is 
much less than the speed of light, with a good approxi- 
mation we can describe electromagnetic interactions be- 
tween charge carriers by instantaneous Coulomb poten- 
tial, and neglect the magnetic fields created by them. 

Finally, we should take into account that graphene is 
placed on a substrate with dielectric permittivity e. A 
physical way to account for the substrate would be to 
modify the lattice action only for plaquettes which are 
inside the medium. Such modification might be advan- 
tageous for studying multi-layered graphene or Casimir 
interactions with graphene sheets. However, since in this 
paper we are interested only in the simplest geometry 
in which the substrate fills half of the three-dimensional 
space, we simply replace e 2 by in all expressions 
above. 

We thus arrive at the following discretized path inte- 
gral representation for the trace Tr em over the states of 
the electromagnetic field in the partition function (10): 



Tr em exp [—H/T^ = J d<t>(s,£,r,z) exp(-S em t, 2)]) 



S,£,T,Z 

Lr-l 

x cxp I -H tb AT + i^2q(s,£_) <f>(s,£,T,z = 0) 

r/Ar=0 \ s,f 



(29) 
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where the lattice action for the electrostatic potential cf> (s, £, r, z) is 



Sem [<j> (S,£,T,Z)] 



(30) 



and 



V3Az e+1 V~3 ( Az\ ( Ko) e+1 



4ire 2 Ar 2 4?re 2 V a / («At) 2 
Sv^a 2 e + 1 3^3 /Az\ _1 («a) e + 1 



167re 2 AzAr 2 16?re 2 \ a / (kAt) 2 

I 



(31) 



For further convenience, we have represented the inverse 
lattice coupling constants fthex and (3 Z in terms of dimen- 
sionless combinations of lattice parameters Az/a, kAt 
and en = 1.946 • 1CT 3 . 

Since the fluctuations of the spatial components 
8b (£, z) and 8 Z (s, £, z) of lattice gauge field can be ne- 
glected, one can also remove the operators exp (±i0xy 
from the tight-binding Hamiltonian (9) in (29). 



> 

CD 




FIG. 4: A comparison of the interaction potential of two 
charges on the 24 x 24 x 24 lattice obtained from the dis- 
cretized action (30) (points) with the potential obtained from 
the solution of the continuum Laplace equation with two point 
sources with opposite charges on the torus of appropriate size 
and with the infinite-space Coulomb potential V (r) = e 2 /r 
(solid lines). The inset shows the level surfaces of electrostatic 
potential on the hexagonal lattice. 

From (30) one can see that the field </> (s, £, z, r) is non- 
propagating, thus it does not describe any dynamics of 
the electromagnetic field, but only the electrostatic in- 
teraction. The interaction potential of two charges at 
distance r is now different from the Coulomb potential 
e 2 /r due to discretization errors of order O (a 2 /r 2 ) and 
O (a Az/r 2 ) . These errors can be systematically reduced 
by constructing improved actions with finite differences 
which involve not only nearest neighbors, but also lattice 



sites separated by two and more lattice spacings. 

On Fig. 4 we compare the interaction potential of two 
static charges obtained from the discretized action (30) 
on the 24 x 24 x 24 lattice with Az = a with the po- 
tential obtained from the solution of the Laplace equa- 
tion with two point sources with opposite charges on the 
torus of appropriate size (-\/3 24a) x (3/2 24 a) x (24 a), as 
well as with the Coulomb potential in the infinite space 

V (r) = e 2 jr. The inset illustrates the level surfaces of 
electrostatic potential on the hexagonal lattice. Fig. 4 
shows that our lattice discretization indeed reproduces 
the continuum electrostatic potential with a good preci- 
sion. For all lattices which we have used for simulations 
discretization errors do not exceed few percents. Our 
lattice regularization of the electrostatic interaction also 
unambiguously fixes the value of the on-site interaction 
potential to u — 26.1 eV. 

According to 1 ' , the bare Coulomb interaction with 

V (r) = e 2 /r between electrons on ir orbitals is in fact 
additionally screened by a factor ~ 2 even in suspended 
graphene due to the influence of electrons on other or- 
bitals. Such screening can be roughly accounted for 
by multiplying the coupling constants (31) by this fac- 
tor. Note that in this case our value of the screened 
on-site interaction potential u' ~ uq/2 is quite close to 
the value u m 10 cV obtained in 4 '. However, a consis- 
tent treatment of this screening of Coulomb potential re- 
quires many technical complications, which are certainly 
beyond the approximations used in this paper. For this 
reason, here we do not take it into account. Wc further 
discuss possible effect of such screening in the concluding 
Section 4. 



1.4. Lattice action for the fermion fields 

We now take the remaining trace over the states of 
the fermionic field in the partition function (10). To this 
end we insert the decomposition of the identity opera- 
tor in the fermionic Hilbcrt space between the operator 
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exponentials in (29): 

1= IJ J dfj a (s,£,T) dri a (s,£,T) 

e -fl.W,r)^(.Ar) \vAs,Z,T)}(r)*(s,Z,T)\. (32) 

Here |t7 ct (s, £, r)) are the fermionic coherent states 42 de- 
fined in terms of the Grassman- valued field r\ a (s, £, r): 

M«, 6 r)> = cxp f £ r? CT (,s, C r) ^ (*, J |0) . (33) 

In the coordinates (s,£) introduced above, the tight- 
binding Hamiltonian (9) can be written as 

k *= E E (34) 

I 



where we have again replaced the abstract site indices 
X, y with the coordinates (s, £) introduced in Subsection 
1.2. h s £; S i(:> are the matrix elements of the single-particle 
Hamiltonian, which acts on single-particle wave functions 
as follows: 

2 

[hi)] {a, £) = -K S jT*l)(p^ + Pb )+mil) (a, £) 

2 

M 08, = -k ]T 1> M - Pb) - m V 08. ■ (35) 

b=0 



We consider now the matrix element of one of the oper- 
ator exponentials in (29) between the states \rj a (s,£,r)) 
and \r\ a (s, £, r + At)) and apply the identity lJ 



(f,| exp \T A.j^cj \n') = exp ( £ (e A ) . . 77, qjj . (36) 
Including also the electrostatic potential 0(s,£,t) introduced above, we obtain 

{v*(s,S,t)\ exp -At # («', O + 

\ cr,s,£,s',£' 

]T ±<^ (5, f , r, z = 0) (s, & (s, K («, £, t + At)) = 
= exp 2 fj a {s,^T)[e X p(-hAT±i ( f>(T))](s,^ S , ,OvAs'^',r + AT)\, (37) 

V" - 1 - c / 

where exp (— hAr ± (t)) is the matrix exponent of the onc-particle operator h(s,£;s',£') At ± 
i<f)(s,£, t, z = 0) 5 (s, £; s',£'). Since the fields Vv(s,£) have opposite charges for er =f and er =4., we take the 
plus sign before the electrostatic potential <j> (s, £,t,z = 0) in (37) for a =f and the minus sign - for a =\.. 
Evaluating this exponential to the first order in At, we obtain: 

[cxp (-hAr ± i<f, (t))] (*, C; s', £') = 6 SS , 6 (£, £') e ±<*(».*.r,*=o) _ 
1 

-At / due ±til - u) ^ s ' i - T - z = 0) h(s,^s\Oe ±lu ^ s '' i '- T ' z= °y (38) 



Here we have used the matrix identity 

where denotes the path-ordering of the second expo- 
nential with respect to the u integration variable. As 



discussed in Appendix B, within the approximations that 
we make in this work one can replace integral over u in 
(38) by a value of the integrand at any u <G [0,1]. For 
definitcness, we approximate the matrix exponential in 
(37) as 

[exp(-/iAT±^(T))] (s, = 
= e ±*C«.^.*=°> (S ss , S (£, C) -Arh (s, £; £')) (40) 
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Note that such a choice is different from the approxima- 
tion discussed in 37 ' 38 . 

Now we insert the approximation (40) into the matrix 
element (37). Finally, we should take the product of such 



matrix elements for all r to obtain the weight of the inte- 
gral over r\ (s, £, r). Taking into account the form of the 
Hamiltonian h (s, £; s' , £'), we obtain the following lattice 
action for the fermion fields 



Stb Via (b,€,t)] = {s,£,T)M a [s,Z,,t;s' ,T']r) a (s',£',t) = 

= ^( S ,e,r)(^( S ,C,r)-e ± ^^ z =°^ (T ( S ,^r + Ar,z = 0)) + 

+ K At J2 vA^ti,T)e ± ^ a ^=V Va ((3,t + p b ,T + AT) + 

+m At i~h {a, £, r) e ±^(«.€.^=<V ( a> £, r + At) - 
-m At ^ ^ 08, t) e ^^=% a (fi, £, t + At) . (41) 

i 



Here we have introduced the fermion hopping matri- 
ces M a with matrix elements M a [s,£, r; s',£',r'], which 
arc the functions of the electrostatic potential field 
<)> (s, £,t,z = 0) in the graphcnc plane. A crucial obser- 
vation is that due to the symmetry between particles and 
holes (which correspond to two different components of 
spin cr with our choice of the ground state (2)) the ma- 
trices Afj- and Afj, are complex conjugate: 

M i (s, t; s', t') = M t (s, £, t; s\ t') (42) 

We note here that in contrast to fermionic actions com- 
monly used in the context of lattice gauge theories, such 
as staggered fcrmions 42,48 , our fermionic action (41) does 
not suffer from the doubling of fermion flavors (see Ap- 
pendix A for the proof). The reason is that we use 
the non-symmetric discretization of the time derivative 
d T V* (s, 6, t) « (T] a (s, £, t + At) - r} a (s, £, t)) /At. For 
lattice discretizations of relativistic field theories, such 



lattice derivative would violate cubic symmetry group of 
the lattice, which is the remainder of the Lorentz invari- 
ance. However, in our case there is no Lorentz invariance, 
and Euclidean time and spatial coordinates enter the ac- 
tion in essentially different ways. Thus we do not break 
any symmetry by using the non-symmetric finite differ- 
ence for the lattice derivative. Interestingly, a similar 
path integral representation of the partition function of 
the tight-binding model (9) with the symmetric lattice 
derivative in the time direction has been considered re- 
cently in 19 , and the two fermionic doublers which appear 
due to such discretization were interpreted as the two 
non-relativistic spin components with cr =t, i- 

Finally, integrating over the fields r\ a (s,£, r) and tak- 
ing into account the relation between fermion hopping 
matrices (42), we arrive at the following representation 
of the partition function (10) in terms of the lattice path 
integral over the electrostatic potential field <f> [s, £, t, z]: 



Z = J Vt)„ (s, £, t) Vr\ a (s, £, t) V<f) (s, £, t, z) exp (- ^ V<rM a [(j> (s, f , t, z = 



Tla - Sem [<j>(s,Z,T,z)]\ = 

V(f> (s, t, z) |det (M t [<f> (s, £,z = Q, r)]) | 2 exp (S em [<f> (s, t, z)}) (43) 

I 

As usual, taking the trace over fermionic states involves rj a (s, £, t), thus anti-periodic boundary conditions in Eu- 
one additional permutation of Grassman- valued fields clidcan time r with period (fcT) -1 should be imposed 
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on them. In practice, this amounts to the replace- 


kAt 




T 


ment <f> (s, £, t, z = 0) — > <fi ( s : £) T i z = 0) + 7r within the 
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fcrmionic part of the action in (43) on a single time slice 
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1.5. Lattice Monte-Carlo simulations 
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Path integral weight in the partition function (43) is 
manifestly positive, thus functional integration can be 
performed numerically by a Monte-Carlo method. Con- 
figurations of the electrostatic potential field 4> (s, £, t, z) 
should be therefore sampled with the weight (43). 
Since this weight includes nonlocal determinant of the 
fermion hopping matrix M a [</)(s,£,t,z = 0)], the most 
suitable simulation method is the Hybrid Monte-Carlo 
algorithm 42,48 . 

We use the so-called ^-algorithm 48,49 , in which 
the squared modulus of the determinant of 
M a [<j> (s, £, t, z = 0)] in (43) is represented in terms 
of the complex- valued pseudo- fermion field \ ( s i T ) : 



|det {M\ | 



VxVx exp 



1 



M t M\ 



(44) 



At the beginning of each Molecular Dynamics trajectory, 
we generate the random pseudo-fcrmion field x accord- 
ing to the weight P[x] ~ exp ( ~x (M^Mf\ x] an d 



then perform the Molecular Dynamics evolution of the 
electrostatic potential field <f> (s, £, t, z) with the force 



F[4>] (s,£,t,z) 








■ s em [4>] 



£z,0 X 



>(s,£,T,Z = 0) 



0(s,£,t,z) 

^M t [4>] M\ [(/>}) ~ L x- (45) 



The corresponding equations of motion are solved by us- 
ing the Sexton- Wcingartcn integrator 48 ' 00 . In order to 
improve the ergodicity of the algorithm, the number of 
integrator steps is drawn from the Poisson distribution 
(see, e.g. ' 1 ). The mean value of this distribution is au- 
tomatically tuned during the thermalization process so 
that the acceptance rate of the algorithm lies in the 
range 0.6 ... 0.9. The integrator step size is then changed 
in such a way that the total trajectory length is equal 
to one. We use the standard Conjugate Gradient algo- 
rithm to invert the operator M-j- [</>] Ml [(/)] in (45). After 
the Molecular Dynamics evolution we perform the usual 
accept-reject step. All random numbers are generated 
by using the ranlux random number generator with 
double precision. 

We also additionally speed up our algorithm by 
applying local heatbath updates 42 ' 48 to the variables 
<fi(s,£,T,z) with z between Hybrid Monte- Carlo 
updates. For both updates, the path integral weight 
(43) is the stationary probability distribution. In ad- 
dition, Hybrid Monte-Carlo updates satisfy the detailed 



TABLE I: The parameters of lattices which we have used in 
our simulations. 



balance condition, and heatbath updates satisfy the lo- 
cal detailed balance 12 . By combining the corresponding 
transition probabilities it is easy to see that the path 
integral weight (43) is still the stationary probability dis- 
tribution for the successive application of both updates, 
despite the fact that the detailed balance condition is 
no longer satisfied. We perform 20 global heatbath up- 
dates between successive Hybrid Monte-Carlo updates. 
Within each global heatbath update, we select at ran- 
dom L x x L y x L T x (L z — 1) lattice sites outside of 
the graphene plane and apply local heatbath updates to 
them. This procedure, while consuming less than 10% 
of the total CPU time, significantly decreases the auto- 
correlation time of the algorithm. We have estimated 
the latter for the physical observables such as the mean 
plaquette, the Polyakov loop and the chiral condensate 
as well as for purely algorithmic parameters such as the 
number of iterations of the CG algorithm and the en- 
ergy difference for the Molecular Dynamics trajectories. 
We have found that for all observables and for all lattice 
parameters which we have used the autocorrelation time 
does not exceed 5 full Monte-Carlo updates (which com- 
prise both Hybrid Monte Carlo and heatbath updates). 



1.6. The choice of lattice parameters 



In practice, the simulations can only be performed for 
the finite lattice sizes L x , L y , L z , L T and at finite nonzero 
values of Ar and to. The results should be then extrap- 
olated to the limits Ar — > with fixed temperature in 
physical units T = (L T Ar) 1 and m — > 0, Az — > 0, 
L x — > oo, L y — > oo, L z — > oo. The latter limit should be 



taken in such a way that (mL x ) 



0, (jnLy 



0. 



(AzL z ) _1 — >■ 0. These limits are completely analogous 
to thermodynamic and chiral limits in lattice QCD sim- 
ulations. 

The parameters of lattices which we have used in 
our simulations are summarized in Table I. In order to 
locate the low-temperature phase of the tight-binding 
model (9), we have considered three different temper- 
atures (T/k = 0.28, T/k = 0.42 and T/k = 0.56). As 
discussed in the Introduction, as long as we consider only 
phenomena which involve electronic degrees of freedom 
and which are characterized by typical energies of or- 
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der of k, the low-tcmperaturc limit can be studied with 
a good precision even by considering the temperatures 
which arc much higher than the room temperature. 

Finite- volume effects are controlled by performing sim- 
ulations at fixed temperature in physical units but at 
different spatial volumes (18 3 and 24 3 ). We use lattices 
with sizes which arc multiples of three, because the Dirac 
points are only covered by discrete lattice momenta on 
such lattices (see Appendix A). For all lattices we also 
assume that Az = a. Indeed, since the Coulomb interac- 
tion potential is anyway approximated with an error of 
order a 2 , it does not make sense to take Az <C a. 

For each set of lattice parameters described in Table 
I, we consider four different values of the staggered po- 
tential m (which also plays the role of the Dirac mass 
at low energies): m/re = 0.1, mjn = 0.2, m/n = 0.3 
and m/n = 0.5. With such values of m, the induced 
gap in the spectrum of the free tight-binding model (9) 
is still much smaller than the energy scale E ~ k at 
which deviations from the low-energy linear dispersion 
relation E (k) = vp k become important (see Fig. 1). On 
the other hand, with such choice of parameters the gap 
width is comparable to the temperature, thus one can ex- 
pect that finite-temperature effects might be quite signif- 
icant. For fixed values of k At and to/k, we change the 
strength of Coulomb interaction by adjusting the cou- 
pling constants (3 Z and fthex in (30) according to (31). 
We consider the values of substrate dielectric permittiv- 
ity uniformly covering the range e = 1.0. ..10.0 with step 
Ae = 0.5. For each data point at fixed values of e, k At, 
L x , Ly, L z , L T and to we have generated 100 statistically 
independent configurations of the electrostatic potential 
field 4> ( s , t, z). 
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FIG. 5: Total CPU time (for a 2.4 GHz Intel Xeon CPU) 
required for one Hybrid Monte-Carlo update as a function 
of the substrate dielectric permittivity e for different lattice 
parameters. 

To illustrate the performance of our algorithm, on 
Fig. 5 we plot the total CPU time required for one Hy- 
brid Monte-Carlo update (plus 20 heatbath updates) as 
a function of substrate dielectric permittivity e in (31) 



for different lattices and for different values of the Dirac 
mass to. One can see that the algorithm significantly 
slows down as we move to smaller e and to, that is, deeper 
in the non-perturbative regime with large coupling con- 
stant and small energy gap. The situation is similar to 
that in lattice QCD, where simulations at small quark 
masses also suffer from significant slow-down sometimes 
called the "Berlin Wall" 53 . 



2. SPONTANEOUS BREAKING OF 
SUBLATTICE SYMMETRY 

2.1. Basic definitions and lattice observables 

In this Section we study the spontaneous breaking of 
sublattice symmetry within the tight-binding model of 
graphene (9) with electromagnetic interactions. Before 
discussing the relevant order parameters, let us consider 
more closely the symmetries of this model. 

In the absence of interactions, the Hamiltonian (9) has 
a global U (2) flavor symmetry. For graphene at half- 
filling, it is explicitly broken down to U {1)®U (1) by the 
Coulomb interaction term. This U (1) (g> U (1) symmetry 
ensures the conservation of the total numbers of charge 
carriers with different spins and cannot be broken neither 
by the staggered potential nor by the Coulomb interac- 
tions. In contrast, sublattice symmetry is a discrete sym- 
metry, which can be realized as reflections with respect 
to the planes which are perpendicular to one of the basis 
vectors e a (see (15)) of the hexagonal lattice and which 
intersect lattice links in the direction a in the middle. It 
can be broken either explicitly, by introducing the stag- 
gered potential, or spontaneously, due to a strong enough 
Coulomb interaction. Close to the Dirac points this dis- 
crete symmetry is enhanced to a continuous chiral sym- 
metry of Dirac fermions, so that the overall global sym- 
metry group is enhanced to U (4). Thus, strictly speak- 
ing, Goldstone's theorem is not applicable to spontaneous 
breaking of sublattice symmetry, but Goldstone bosons 
might still appear as effective degrees of freedom at low 
energies . Finally, we note that the symmetries of the 
tight-binding model arc quite different from those of the 
staggered fermionic action, which was used for numeri- 
cal simulations of the effective field theory of graphene 
in . Staggered fermions have a U (1) symmetry asso- 
ciated with charge conservation as well as a U (1) chiral 
symmetry which is explicitly broken by the mass term, 
thus for staggered fermions discrete sublattice symmetry 
is replaced by a continuous symmetry at all energy scales. 
On the other hand, only the total charge of both flavors 
is conserved, but not the charges of each flavor 29,35 ' 40,41 . 

An obvious order parameter for the spontaneous break- 
ing of the discrete sublattice symmetry is the difference 
of the particle number densities on the two sublattices of 
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the hexagonal lattice A at: 



Aft 



Ajv) = - 



ZL X Ly 



T dlogZ 
L x L y dm 

Tr fAve-^ 



A ^ = E («' ^ - # 3 , & 0) (46) 



'A 



^Re(Af t - 1 (a^,r;a,e,r))- 



Z/x Lj, i r 



^Re(A/f 1 (/3,e,T;/3,C,r)),(48) 



A simple calculation within the Dirac approximation 
shows that in the absence of interactions and at suffi- 
ciently small values of the staggered potential m ( A at ) 
is a linear function of m. It is also convenient to introduce 
the susceptibility of A at as 



d(A 



Xn 



dm 



m-K) 



(47) 



At sufficiently small temperatures, when only the linear 
part of the spectrum contributes to the expectation val- 
ues, ( A ) jv and xn can be expressed in terms of the chiral 
condensates and chiral susceptibilities of the two flavors 
of Dirac quasiparticles. By analogy with chiral symme- 
try breaking in gauge theories, one can expect that at 
a second-order phase transition the susceptibility (47) 
should diverge. 

To obtain the expression for the expectation 
value (46) on the lattice, we insert the operator 

Ajv = E te i>a (a, - fa (A {fi, 0) be- 

tween O'th and L T 'th factors in the Feynman-Kac rep- 
resentation (29) of the partition function (10). After 
bringing the operators t\i a (s,£) and ip^ (s, £) to the nor- 
mal order, the integral over the fermionic coherent states 
can be easily taken. We then obtain the following expres- 
sion for ( Ajy ) in terms of the fermionic hopping matrix 



where the brackets (...) on the r.h.s. denote averaging 
over the electrostatic potential field <j> (s, t, z) with the 
weight (43) and AfT 1 is understood as the matrix inver- 
sion of the fermion hopping matrix Af-j- in (41). 

By analogy with chiral condensate measurements in 
lattice QCD, one can expect that the expectation value 
( Ajv ) in the limit m — s- should be distorted by large 
finite-volume corrections as well as by systematic errors 
of extrapolation torn = 0. The reason is that the "chiral" 
limit m — > and the thermodynamic limit L x — > oo, 
L y — >• oo do not commute (see e.g. 48 , Sec. 15.2.1), and 
strictly speaking ( A N ) at m = should be zero in any 
finite volume. The numerical value of the condensate 
therefore strongly depends on the way of extrapolating 
it to the limit m — > 0. Therefore we also consider the 
dispersion of A at, which is finite in finite volume and is 
thus much less affected by these numerical uncertainties: 



«A 2 ^» 



1 



2 L X Ly 



-Tr 



(A^e-^) 



L^Lii ( A 



(49) 



Similarly to the susceptibility xn- (( A^ )) should signifi- 
cantly increase or diverge at the phase transition. The ex- 
pression (49) contains contributions both from connected 
and from disconnected fermionic diagrams: 



((A^)) = ((A^)) conn . + ((A^)) di 

(( a^ >u. = - L 2 - Re ( E < M f 1 ( s > r > s > t ) ) 

X V T £,T s 



( A/f 1 (s, t r; s, r) Mf 1 {a, r; s, r) } M f' («> ^ ^ P, M t* (/?> T ' Q ' ^ T ) i 



■2(Mjr 1 G8,e,r;a,^,r) Mf 1 (a, r; (3, £, r) ) ) , 



(( a 2 n )) dlsc . = - - - E M E ( M f 1 («> r ; a ' r ) - M f 1 OS- r ; ^ T ) 

x y t r 



4 



L x L y L T 



- E ( M f 1 ( a > f » r ; a ' s> r ) - M f 1 ' & r ; & & r )) > 2 > ( 5 °) 



where again the brackets (...) on the r.h.s. denote aver- aging over the electrostatic potential field with the weight 
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(43) and the first summand in (( Aj^ }} conn , arises due to 
an additional exchange of fcrmion operators ip a (s, £) and 
l ijr a (s, £) which is necessary in order to bring them to the 
normal order within the expectation value of the four- 
fermion operator. We have found that the disconnected 
contribution is much noisier than the connected one (see 
Fig. 10 below), so that the computer time required to es- 
timate it with sufficient precision becomes prohibitively 
large. For this reason, we have considered only the con- 
nected part ((A N )) C onn.- It should be also noticed that 
the connected contribution (( )) coran , does not in gen- 
eral correspond to the expectation value of any opera- 
tor in canonical formalism, that is, it cannot be repre- 
sented in the form Z _1 Tr (Oe^ 13 ^ with some opera- 
tor O. This is in contrast with lattice QCD, where the 
connected contributions to the correlators of two-fermion 
operators can be interpreted in terms of physical meson 
states consisting of quarks with different flavors. The 
reason is that in our case the two fermion flavours have 
opposite charges and thus couple differently to the elec- 
trostatic potential field (f>. 



2.2. Simulation results 

On Fig. 6 we plot the expectation value A at of the 
difference of particle number densities on even and odd 
lattice sites as a function of the substrate dielectric per- 
mittivity e at fixed values of mj K (plots on the left) and 
as a function of m/ K at fixed values of e. We present the 
results for the 24 4 lattice at the temperature T/k = 0.56 
(kAt = 0.075, plots at the top) and at T/k = 0.28 
(k At = 0.15, plots at the bottom). (An) gradually in- 
creases as we move into the strong-coupling region (small 
e) or to larger values of mj K. In the weak-coupling region 
(large e) ( Ajv ) is almost a linear function of 771, while in 
the strong-coupling limit this dependence becomes essen- 
tially nonlinear. 

In order to extrapolate (An) to the limit m — > 0, we 
fit the dependence of (An) on m/n at fixed e with a 
quadratic polynomial and use the value of this polyno- 
mial aim — as an estimate of ( Ajv ) |m-K>- The result of 
such extrapolation and the corresponding fits are shown 
on Fig. 6 with solid lines. All fits have /d.o.f. of or- 
der unity. On Fig. 7 we also compare the dependence of 
the extrapolated values of ( An ) on e for different lattice 
parameters. While at higher temperature (T/k = 0.56) 
the result of extrapolation is equal to zero within error 
range, at T/k = 0.42 and T/k = 0.28 one can clearly see 
that ( Ajv ) remains finite in the limit m — > for e < 4 
and grows as the temperature decreases. This indicates 
that sublattice symmetry of the tight-binding model (9) 
is spontaneously broken due to Coulomb interaction at 
T < 0.42 k and at e < 4. For the 24 4 lattice at T/k = 0.28 
the extrapolated value An is somewhat larger than for 
the 18 4 lattice at the same temperature, as it should be 
for spontaneous symmetry breaking. 



In order to get further insight into the nature of the 
transition to the spontaneously broken phase, on Fig. 
8 we plot the susceptibility xn as a function of sub- 
strate dielectric permittivity e for different lattice pa- 
rameters. The susceptibility was also obtained from the 
quadratic fits of the dependence oi (An) on mj k as the 
first derivative of the fitting polynomial at m — 0. At 
T/k = 0.56 and T/k = 0.42 xn monotonically grows as 
e decreases, reaching its maximal value at e = 1, that 
is, for the strongest Coulomb interaction. In contrast, at 
T/k = 0.28 xn becomes a non-monotonic function of e 
with a characteristic peak at e w 4. For the 24 4 lattice 
this peak is somewhat sharper than for the 18 4 lattice. 

To present an additional evidence of the existence of 
this peak which is independent of any fitting procedure, 
on Fig. 9 we plot the connected part of the dispersion 
of the difference of particle number densities on two sim- 
ple sublattices ((A N )) conn which was directly calcu- 
lated on different lattices according to (50) at m = 0.2 k. 
((A N )) conn , as a function of e also has a distinct peak 
at 4 < e < 5 for T/k = 0.28, and a somewhat less pro- 
nounced peak at 3 < e < 4 at T/k = 0.42. Interestingly, 
for ((A N ))conn. the height of the peaks practically does 
not depend on the lattice size. 

Our reason for considering only the connected part of 
((A N )) is that the disconnected part turns out to be 
much noisier. In order to illustrate this observation, on 
Fig. 10 we plot both the connected and disconnected 
contributions for the 18 4 lattice with kAt = 0.2 and 
m/re = 0.2, 0.3, 0.5 as a function of e. The disconnected 
contribution was calculated using 500 Gaussian stochas- 
tic estimators (see e.g. 48 , Sec. 11.1), which required sev- 
eral hundreds core-hours per data point (for a 2.4 GHz 
Intel Xeon CPU) in the strong-coupling regime and at 
the smalles value of m (m/ k = 0.1). While this amount 
of computer time is already close to the time required 
to generate the field configurations, the numerical errors 
of the disconnected part are still much larger than those 
of the connected part, especially in the strong-coupling 
regime. Thus the reliable estimate of the former would be 
prohibitively expensive, and we do not consider it here. 
Probably some progress in this direction could be made 
by using graphic cards (GPUs) or by applying some 
more refined measurement procedure. We leave such de- 
velopments as a direction for further work. Since we only 
use (( A^r )) for a more precise location of the transition 
point, one can expect that the connected contribution 
alone can also be a good estimator. Indeed, if there is a 
second-order phase transition in the tight-binding model 
(9), than all order parameters have singularities at a sin- 
gle value of e, and if there is a crossover, than the critical 
value is not well-defined anyway. 

We conclude that the behavior of Ajv and xn at 
T = 0.28 k is suggestive of a second-order quantum phase 
transition with respect to substrate dielectric permittiv- 
ity e at the critical value 4 < e c < 5. At T = 0.56 k sub- 
lattice symmetry is always restored in the limit to — > 
and this phase transition is obviously absent. The inter- 
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FIG. 6: Differences of particle number densities on simple sublattices in graphene as a function of substrate dielectric permit- 
tivity e (on the left) and staggered potential m (on the right) on the 24 4 lattice. Above: at T/k = 0.56 (k At = 0.075). Below: 
at T/k = 0.28 (k At = 0.15). Points with solid lines through them on the plots on the left are the results of extrapolation to 
the limit m — s> 0. These solid lines are weighted splines and are shown to guide the eye. Solid lines on the plots on the right 
are the quadratic fits which were used for extrapolation. 



pretation of the lattice data at T = 0.42 k is not quite 
straightforward: while the extrapolation of An to the 
limit m — > yields nonzero result and the dispersion 
(( A N )) CO nn. as a function of e still has the characteristic 
peak, the susceptibility xn is a monotonic function of e. 
It is likely that at this temperature we arc close to the 
cndpoint of the second-order phase transition line in the 
parametric space (e,T), where the second-order phase 
transition disappears or turns into a crossover. The fact 
that at T = 0.42k the peak of (( A N )) conn , is situated 
at smaller e than at T = 0.28 k suggests that for higher 
temperatures the critical value of e becomes somewhat 
smaller. As discussed in the Introduction, in this pa- 
per we are interested in the low-temperature phase of 
the theory. Therefore we conclude that the expected 
pattern of spontaneous symmetry breaking in the low- 
temperature phase is observed for T < 0.28 k, and do not 
study the presumable finite-temperature phase transition 
at T = T c ps 0.42 n. Finally, we note that the volume- 
independence of (( Ajy- }} co7m . at m/n = 0.1 might indi- 
cate that for nonzero m the second-order phase transition 



at T < 0.28 k and e w 4 also turns into a crossover. 



Combining all the data, we estimate the critical value 
of the substrate dielectric permittivity as e c — 4± 1. Ex- 
pressing the effective QED coupling constant a in terms 
of e as a — — -^r, an w 1/137, we find that this 
value corresponds to a c = 0.9 ± 0.2. This estimate is 
in agreement with the results of simulations of graphene 
effective field theory with staggered fcrmions J<, ~ 28 ' 30,34 , 
where the second-order phase transition to the phase with 
spontaneously broken chiral symmetry was observed at 
a c = 1.11 ± 0.06. We further discuss the phase structure 
of the tight-binding model (9) in the concluding Section 
4. 
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FIG. 7: Extrapolation of the difference of particle number 
densities on simple sublattices Ajv to the limit m — > as 
a function of substrate dielectric permittivity e at different 
lattice parameters. Solid lines are the weighted splines which 
are plotted to guide the eye. 
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FIG. 9: Connected part of the dispersion of the differ- 
ence of particle number densities on two simple sublattices 
({ A% )} C onn. as a function of substrate dielectric permittiv- 
ity e at m/n = 0.1 on different lattices. Solid lines are the 
weighted splines which are plotted to guide the eye. 




FIG. 8: Susceptibility xn of the difference of particle number 
densities on simple sublattices Ajv as a function of substrate 
dielectric permittivity e at different lattice parameters. Solid 
lines are the weighted splines which are plotted to guide the 
eye. 



3. GRAPHENE CONDUCTIVITY FROM THE 
GREEN-KUBO RELATIONS 



3.1. Basic definitions and lattice observables 



In this Section we study numerically the conductivity 
of graphene monolayer, that is, the linear response of the 
electric current to the applied homogeneous electric field. 
In order to define the operator of electric current within 
the tight-binding model (9), we consider the time evolu- 
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FIG. 10: A comparison of connected and disconnected contri- 
butions to the dispersion of the difference of particle number 
densities (50) on simple sublattices ({ }) for the 18 4 lattice 
at ftAr = 0.2. 



the form 



(51) 



(4). This leads to the charge conservation equation of 



where t is the real (Minkowski) time and J a (s,£) is the 
operator of the electric current flowing through the lattice 
link which goes in direction a and originates from lattice 
site with coordinates (s, £). It is equal to the difference 
of the currents J a , a (s,£) of "particles" and "holes": 

Ja.b (0 = i>l (J3, £ + Pb ) e=F^«)^ ( a , £) - 

j b (a, ee J b (0 , J b (0, ee J b (e - Pb ) . (52) 
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To study the conductivity of graphene within the linear 
response theory, we have to introduce the classical time- 
dependent background electromagnetic field in the tight- 
binding Hamiltonian (9). This amounts to replacing the 
operators 9xy in (34) with the corresponding classical 
variables (16). As it should be, the electric current oper- 
ator J b (s, £) is then equal to the derivative of the Hamil- 
tonian (9) over the classical link variables 9 b (s,£). Cor- 
respondingly, the linear response of the electric current 
to a small variation 59 b (£, t) of the classical link variables 
is given by 

(J 6 (£,*))= J] j dt'G Rbc (ZM'J) SOc (£',*')• (53) 

Here J b (£, t) is the current operator in the Heisen- 
berg representation and G R bc (£ 5 1; t') is the retarded 
current-current correlator 



xTr 



,-H/T 



(54) 



with 9 (t — t') denoting the Heaviside step function. 

Let us now consider the infinitesimal spatially homo- 
geneous time-dependent electric field 8E (t) = 5 A (t) . 
According to the definition (16), the corresponding 
variation of the classical link variables is 59b (£, t) = 

a (e b -5A(t)^j. Performing the Fourier transform 

6A(w) = 5E(w)/w = J dte- lwt 5A(t) and taking into 
account the spatial homogeneity, we can also write the 
relation (53) as 



J b (w) > = - V G Rbc (w) U ■ SE (w) 



(55) 



where 



G Rb c (w) 



E 

6 . 



dt( 



G Rbc (0,0;U)- (56) 



The conductivity of graphene is defined as the coeffi- 
cient relating the total charge transported through unit 
length per unit time and the applied electric field 1 ' 9 ' 10 ' 55 . 
Since the canonical dimensionality of the electric field 
strength is L~ 2 (where L is the unit length), a (w) is a 
dimensionless quantity. For conversion to the SI system 
of units, it should be multiplied by e 2 /ft. 

For simplicity we assume that the electric field 5E (w) 
is parallel to one of the lattice link vectors e bo . Taking 
into account that the side of the plaquette of the dual lat- 
tice which is perpendicular to e bo is equal to y/3 a and av- 
eraging over all equivalent directions bo, we arrive at the 
following expression for the AC conductivity of graphene: 



a (w) 



G R bc (w) T bc 
3 a/3 w 



(57) 



where 



T, 



bc 



e b -e c = 3/2 6 bc -l/2 



(58) 



and we assume summation over the repeated indices b, c. 

In practice, the AC conductivity a (w) can be ex- 
tracted from the Euclidean current-current correlator, 
which we define as 



G(r) 



3v^3 L x L y 



,-H/T 



x Tr [e™ J b (0 e-™ J c (£') e 
with the help of the Green-Kubo relations 56-59 : 

oc 

dw 



G(r) 



2tt 



K (w, t) a (w) 



where the thermal kernel K (w, r) is 59 

2w cosh (w (r — ^)) 



K (w, t) 



sinh(#) 



(59) 



(60) 



(61) 



In Appendix C we derive explicit expressions for the 
Euclidean current-current correlator G^ (r) and the AC 
conductivity (w) for the tight-binding model (9) in 
the absence of Coulomb interaction. For w -C k, a^ 0>> (w) 
can be approximated as 



(T (0 » (tl))«Si(tl)) + 



+ ■ 



1 (w — 2m) ( 4m 



tanh ( — ) , (62) 



with 3 being some constant. The 5- function singular- 
ity at w = is a common feature of all ideal crystals 
which arises due to the absence of scattering of charge 
carriers. Thus, strictly speaking, cr^ (w) has no well- 
defined zero-frequency limit. A commonly quoted uni- 
versal value o"o = 1/4 (in units of e 2 /h) is obtained from 
(62) at w > T, w > m (but still w < k) 10 ' 55 . Fre- 
quency dependence of cr^ (w) is illustrated on Fig. 16 
(see Appendix C). 

In order to obtain the expression for the discretized 
correlator (59) on the lattice, we insert the current oper- 
ators (52) between O'th and L T 'th and between (r/Ar)'th 
and (r/Ar + l)'th factors in the Feynman-Kac represen- 
tation (29) of the partition function (10). For the time 
being we assume that t/0. Repeating the derivation of 
the fermionic lattice action presented in Subsection 1.4, 
we arrive at the following expression for the discretized 
correlator (59) in terms of the fermionic path integral: 



G(r) = Z- 1 j 



Vf) a Vri a V<t> ■ 



Tbc 



3\^3L x Ly 

E (°)^,6»7<7 (0)) (f)*' (O-V.cV (r)) x 



x exp 



r) a M a [4>] r)„ - S e 



(63) 
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where we have omitted the arguments of the field vari- 
ables for the sake of brevity and ja.b is the one-particle 
operator of the total current of particles with spin a on 
the whole lattice, which is defined by the identity 

£ 

I 



Correspondingly, the symbol r\ a (r) in (63) denotes the 
one-particle wave function r\ a (s, £, r). 



Integrating over the fermion fields and taking into ac- 
count that Mj, = Mf we obtain the following expression 
for the correlator G (r) in terms of the fermionic hopping 
matrix Mf (s, £, r; s', t'): 



G (r) = ( ^ Tr (j t , 6 Mf 1 (0, r) j t , c Mf 1 (r, 0)) ) + 

+ ( t> ReTr (j t , 6 Mf 1 (0,0)) ReTr (j t , c Mf 1 (r,r)) ), (65) 



where (...) denotes averaging over the electrostatic 
potential field (j>(s,£,T,z) with the weight (43) and 
Mf 1 (t, t') is treated as a one-particle operator. Cor- 
respondingly, the trace in (65) is taken over one-particle 
states. The first and the second summand in (65) are 
the contributions of connected and disconnected fermion 
diagrams, respectively. 

When t = 0, an additional interchange of field opera- 
tors tp a (s, £) and (s' , £') is required in order to bring 
them in the normal order and to apply the Feynman-Kac 
transformation. This leads to an additional contact term 
at r = 0, so that 



G(0) 



2 Tfc c 



(ReTr (j fjb j t , c Mf 1 (0, 0)) ) 



2 Tbc 



ReTr ( j t , 6 Mf 1 (0,0) j t , c Mf 1 (0,0) 1 )(66) 



3.2. Simulation results 
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FIG. 11: A comparison of connected and disconnected con- 
tributions to the Euclidean current-current correlator (59) 
(first and second summands in (65)) for the 18 4 lattice with 
K At = 0.2 and m/n = 0.2. The inset shows both contribu- 
tions to the correlators G (r) for r close to /3/2 in a larger 
scale. 



First we estimate the connected and disconnected parts 
of the correlator (59), that is, the first and the sec- 
ond summands in (65). Both contributions are shown 
on Fig. 11 for the 18 4 lattice with kAt = 0.2 and 
to/ ' k = 0.2. Disconnected contributions were estimated 
using 500 Gaussian stochastic estimators . One can 
readily see that the disconnected contribution is much 
smaller than the connected one, and the relative statis- 
tical errors are much larger. From Fig. 11 one can also 
see that the relative importance of this disconnected con- 
tribution is somewhat higher for r close to /3/2 and/or 
for smaller values of e. As discussed in Subsection 2.2, 
estimating the disconnected contributions with sufficient 
precision in the strong-coupling regime by using our cur- 
rent measurement methods would require prohibitively 
large computer time. For these reasons, we disregard 
them in what follows and leave their detailed study as a 



direction for future investigations. 

Connected contributions to Euclidean current-current 
correlators (59) on the 24 4 lattice at the temperature T = 
0.56 k (kAt = 0.075), m/n = 0.1, 0.5 and T = 0.28k 
(k At = 0.15), to/k = 0.1, 0.5 are plotted on Fig. 12 for 
different values of the substrate dielectric permittivity 
e. One can see that as e decreases and the Coulomb 
interaction becomes stronger, the correlators decay much 
faster, which, according to (60), indicates that the AC 
conductivity a (w) becomes smaller in the low-frequency 
region. This effect becomes more prominent at lower 
temperature or at larger values of the staggered potential 

TO. 

For reference, on Fig. 12 we also plot the current- 
current correlators for the non-interacting tight-binding 
model (sec Appendix C for an explicit expression). The 
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FIG. 12: Euclidean current-current correlators (59) on the 24 4 lattice at different values of substrate dielectric permittivity e. 
Above on the left: for T/k = 0.56 (kAt = 0.075) and m/n = 0.1. Above on the right: for T/k = 0.56 (kAt = 0.075) and 
m/n = 0.5. Below on the left: for T/k = 0.28 (kAt = 0.15) and m/n = 0.1. Below on the right: for T/k = 0.28 (kAt = 0.15) 
and m/n = 0.5. 



free correlator obtained from the expression (C9) with 
continuous Euclidean time r is plotted with black solid 
line. The corresponding numerical result which was cal- 
culated according to (65) and (66) with Mj- given by (41) 
and with </> (s, £, t, z) = is shown with black circles. A 
comparison of the results of these two calculations sug- 
gests that the effect of discretization of Euclidean time r 
on the current-current correlators should be rather small. 

A commonly used method to invert the integral equa- 
tion (60) and to estimate the AC conductivity a (w) from 
the values of G (t) in a discrete set of lattice points is 
the Maximum Entropy Method (MEM) 58 ' 59 . However, 
in practice we have found that for our data MEM does 
not produce stable results for a (w) in the low-frequency 
limit (w < T). In particular, it docs not reproduce the 
free AC conductivity crW (w) when supplied with the free 
Euclidean correlator G^ (r). This fact can be probably 
explained by the singularity and discontinuity of (w) 
at small frequencies, which cannot be reproduced by the 
smooth basis functions used in MEM 58 ' 59 (we have used 
the modified thermal kernel (61) introduced in 59 ). In 
fact, MEM tends to simply smear the function a (w) at 



low frequencies, so that the numerically obtained func- 
tion a (w) is smooth and shows no signatures of the gap. 
In this situation the values of a (w) at low frequencies are 
quite meaningless and cannot be compared to the univer- 
sal limiting value cto = 1/4 10,55 . Thus we conclude that 
MEM does not give a reliable estimate of the AC con- 
ductivity of graphene at small frequencies. The situation 
could be probably improved by more advanced modifica- 
tions and tuning of the method, which are out of scope 
of the present paper. 

In order to obtain an estimate of the conductivity 
which is free of the ambiguities introduced by MEM, let 
us consider the Euclidean correlator (59) at r = P/2. 
According to (63) and (61), its value can be represented 
as 

oo 

The weight factor — ^ — r is finite at w — > and decays 

6 sinh(^) J 

exponentially at w > T. The integral in (67) is thus 
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FIG. 13: Smeared low-frequency conductivity (68) in units of e 2 /h as a function of substrate dielectric permittivity e at different 
values of the ratio m/re. On the left: on 24 4 lattice at T = 0.56 k (kAt = 0.075), on the right: on 24 4 lattice at T = 0.28 k 
(k At = 0.15). Points with solid line through them is the extrapolation to the limit m — ► 0. Solid lines are the weighted splines 
which are plotted to guide the eye. 



saturated in the region with w < T. It is thus natu- 
ral to introduce the conductivity a smeared over small 
frequencies as 



(T =JV- 1 ^ 



dw 2w 



2ir sinh I 



- a (w) 



1 



2T i 



irT 2 



G(/3/2), (68) 



where J\f is the normalization factor: 



TV = 



/ , lJ w t = ttT 2 . Analytical calculation of a within 

^ Sln H2T,) 

the Dirac approximation to the non-interacting tight- 
binding model (see Appendix C) shows that it is quite 
close to the limiting value Uq = 1/4 and does not depend 
on temperature in the limit m — > 0. 

Smeared low-frequency conductivity a for the 24 lat- 
tice at the temperature T = 0.56k (kAt = 0.075) and 
at T = 0.28k (kAt = 0.15) and with different values 
of m/ft is plotted on Fig. 13 as a function of substrate 
dielectric permittivity e. At nonzero to d gradually de- 
creases with e. As to becomes smaller, a becomes al- 
most a constant function for e > 4, but changes faster 
at e < 4. In order to extrapolate the conductivity to the 
limit to — y 0, we fit its dependence on m at fixed e with 
quadratic polynomial and use the value of this polyno- 
mial at m = 0. All the fits yield \ 2 /d.o.f. of order of 
unity. Results of such extrapolation for different lattices 
are summarized on Fig. 14. Solid horizontal lines on 
the plot correspond to the results of analytic calculation 
of in the non-interacting tight-binding model with 
to = (see Appendix C). 

The extrapolated conductivity a is practically constant 
and close to its value in the non-interacting tight-binding 
model at e > 4 for T = 0.56 k and for T = 0.42 k, 
but quickly decreases with e at e < 4. At T = 0.28 k, 
the behavior of the conductivity is essentially the same, 
but the critical value of e at which a starts decreasing 
is somewhat higher, e w 5. One can also note a slight 




FIG. 14: Smeared low-frequency conductivity a in units of 
e 2 /h after extrapolation to the limit m — > as a function 
of substrate dielectric permittivity e for different lattice pa- 
rameters. Solid horizontal lines on the plot correspond to the 
results of analytic calculation of <r' ' in the non-interacting 
tight-binding model with m = 0. Solid lines which are plot- 
ted through data points are the weighted splines which are 
shown to guide the eye. 



decrease in the conductivity at larger e, which is in agree- 
ment with perturbative calculations'' . Remarkably, the 
critical values of e which separate the regimes of con- 
stant and decreasing conductivity coincide with the crit- 
ical values which were obtained in Section 2 from the 
analysis of spontaneous breaking of sublattice symmetry 
at the corresponding temperatures. At e = 1, when the 
strength of Coulomb interaction is maximal, the smeared 
conductivity is still finite and comprises 20 — 30% of its 
weak-coupling value. Such behaviour of the conductivity 
should be contrasted with the results obtained from sim- 
ulations of the graphene lattice effective field theory 34 , 
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for which the low-frequency conductivity in the strong- 
coupling phase decreased by at least one order of mag- 
nitude as compared to its value in the weak-coupling 
phase. It seems therefore that the semimetal-insulator 
phase transition associated with spontaneous symmetry 
breaking is somewhat softer for the tight-binding model 
than for the graphene effective field theory. 

The fact that the numerically obtained weak-coupling 
limit of a (m — > 0) for kAt = 0.15 and ft At = 0.20 de- 
viates from the result of analytical calculation (solid hor- 
izontal line on Fig. 14) suggests that our extrapolation 
procedure tends to overestimate the value of the conduc- 
tivity, with extrapolation error being as large as 10% at 
« At > 0.15. On the other hand, at kAt < 0.1333 nu- 
merical results in the weak-coupling regime agree nicely 
with the analytical result obtained in the non-interacting 
model. 

We also note that in the weak-coupling limit the ex- 
trapolated conductivities a (m — > 0) differ quite signifi- 
cantly for different temperatures, in contrast to the re- 
sult (C12) which was obtained in Appendix C in the 
Dirac approximation. This indicates that the deviations 
from the linear dispersion law E = vf k at E > ft and 
the finite width of the valence band of graphene (with 
Emax = V9ft 2 + m 2 ) are still important at our values of 
the temperature. 



4. DISCUSSION AND CONCLUSIONS 

In this paper we have presented the results of numer- 
ical studies of the tight-binding model of graphene with 
Coulomb interaction. We have assumed that the strength 
of Coulomb interaction is controlled by substrate dielec- 
tric permittivity e, so that the QED coupling constant 
«o ~ 1/137 is multiplied by the factor 2/ (e + 1). 

Our results indicate that for sufficiently strong 
Coulomb interaction, that is, at substrate dielectric per- 
mittivities e < 4, and at sufficiently small temperatures 
(T < 0.28 ft) the symmetry between the two simple sub- 
lattices of the hexagonal lattice of graphene is sponta- 
neously broken by a nonzero expectation value (Ajv) of 
the difference of the numbers of particles localized on the 
sites of each sublattice. At the critical value e « 4, the 
susceptibility \n of ( Ajv ) as well as the connected part 
of the dispersion of Ajy (( A 2 ^ )) conn . have distinct peaks 
indicative of a second-order phase transition. This result 
agrees with the results of simulations of graphene effec- 
tive field theory with staggered Dirac fermions '. Our 
estimate of the critical value e c = 4±1 corresponds to the 
effective QED coupling constant a = ^- = 0.9 ±0.2, 
which agrees with the value a c = 1.11 ± 0.06 obtained 
in 26 ~ 28 ' 30,34 . This fact suggests that lattice artifacts of 
staggered fermions have no significant effect on the posi- 
tion of the semimetal-insulator phase transition. 

At higher temperature (T = 0.42 k) sublattice sym- 
metry is still broken, but the transition to the broken 
phase becomes softer and looks more like a crossover. 



In particular, in this case only ((A 2 v )) colm . has a char- 
acteristic peak, and xn is a monotonic function of e. 
At even higher temperature, T = 0.56 ft, there are no 
signatures of spontaneous symmetry breaking. It seems 
therefore that at T = 0.42 k we are in the vicinity of the 
finite-temperature phase transition at which sublattice 
symmetry is restored. 

In order to quantify the electronic transport proper- 
ties of graphene, we have considered the low-frequency 
conductivity a smeared over frequencies w < T accord- 
ing to (68). We have found that for all temperatures 
which we have considered the conductivity a quickly de- 
creases with e at £ < 4 down to a finite value which 
is around 20 — 30% of its weak-coupling limit for the 
strongest coupling (e = 1, which corresponds to sus- 
pended graphene). It turns out that for the tight-binding 
model the low-frequency conductivity decreases some- 
what slower than for the graphene effective field theory, 
which suggests that the insulator-semiconductor phase 
transition is somewhat weaker in our case. At e > e c , a 
practically does not depend on e and gradually decreases 
with temperature. Such behavior of the conductivity in- 
dicates that the phase transition associated with spon- 
taneous symmetry breaking might persist at higher tem- 
peratures, but become weaker. For example, there could 
be a second-order phase transition at small temperatures 
(T < 0.28 k) and a crossover at higher temperatures. A 
more accurate analysis of finite-temperature and finite- 
volume effects is required in order to classify the order 
of the observed phase transition and to obtain the corre- 
sponding critical exponents. 

Finally, let us comment on the obvious discrepancy be- 
tween the results of lattice simulations (reported in the 
works 26-30 ' 34,35 and in this paper), which suggest that 
suspended graphene should be an insulator, and most re- 
cent experimental results 21 ' 22 , which find no signature of 
an insulating state in suspended graphene. As discussed 
in ' , Coulomb interactions in graphene can be addi- 
tionally screened both by valence electrons and by elec- 
trons on other orbitals of carbon atoms. This screening 
becomes stronger at small momenta and can effectively 
decrease the electromagnetic coupling constant by a fac- 
tor of two or larger. Due to the screening, the critical 
value of the substrate dielectric permittivity might even- 
tually become less than e c = 1, which would make the 
insulator-semimetal phase transition inobservable in the 
real world. Since our simulations automatically take into 
account the screening of Coulomb interactions by valence 
electrons, this effect might be explained by the influence 
of electrons on other orbitals. As well, the size of the lat- 
tices which we use for our simulations might be too small, 
so that the momentum scale at which the screening be- 
comes sufficiently strong is not yet reached' 1 . Another 
possible cause of this discrepancy is the wrong value of 
the on-site interaction potential Uo, which is also a free 
parameter of the tight-binding model. If we couple the 
tight-binding model to the lattice gauge field, as in (9), 
the value of uq is fixed by the choice of the lattice action 
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and might be quite different from its physical value. All 
these conjectures require separate investigations. 
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Appendix A: Eigenspectrum of the tight-binding 
Hamiltonian on the finite lattice 



In this Appendix we discuss the spectra of the free 
single-particle Hamiltonian (35) and the fermion hopping 
matrices M a introduced in (41) on lattices of finite size. 

Invariance of the one-particle Hamiltonian (35) under 
translations implies that its eigenfunctions take the form 



1>C (", & q) = N a £ (?) cxp (i qg), 

A {P, g, q) = A/>, c (g) exp (i qg), 



(Al) 



where we have introduced an additional label £ to distin- 
guish between states with equal momenta but different 
energies. The components of the wave vector q in Carte- 
sian coordinates are 



fc _ gi k _^_qi_ 

V3a V 3a 3a 



(A2) 



It is easy to check that the functions (Al) are the eigen- 
functions of the one-particle Hamiltonian (35) with the 
eigenvalues 



E C (q) = (E (q) = CVm 2 + n 2 [1 (q) | 2 , (A3) 
where ( takes values ( = ±1 and 

$ (q) = J2 e tqPb = 1 + e- iqi+iqz + e^ 91 . (A4) 

b 

The ratio of the normalization coefficients Af a ,C (?) an( i 

■A/>,c (q) is: 



MpX («)/^a lC (?) 



m - E c (q) 



(A5) 



This equation and the normalization condition 
||V(g)|| 2 = L x L y (\jV aX (q)\ 2 + \jV PX (q)\ 2 ) fix the 
values of the normalization coefficients: 



•A/>,c (q) = ~(e 



-i0(q) 



1 E(g) + (m 
2E(q) L x L y ' 

I E(q)-(m 
2E{q) L x L y 



(A6) 



where $ (g) = arg $ (q). 

Dirac points correspond to lattice momenta q with 
$ (q) = 0. This condition is equivalent to the two equa- 
tions 

cosgi + cosg 2 + 1 = 0, sin q\ = — sin q 2 . (A7) 

Solving these equations, we find two Dirac points q^ 

Tj, Linear expansion of the 



(±) 



q^ 



with q\ 

dispersion relation (A3) with m = near these points 
leads to the well-known result E (fc) = vp\k\ with the 
Fermi velocity vp = 3/2 n a. The largest eigenvalue of 
the single-particle Hamiltonian (35) is 



E„ 



E+ (0) = \/m 2 + 9k 2 



(A8) 



We should also take into account the boundary condi- 
tions (13), which constrain possible values of q: 

q\L x = 2nmi, mi G Z 
q2L y - qiLy/2 — 2nm,2, ra 2 eZ. (A9) 

Expressing qi , q 2 in terms of mi and mi we find 
27rmi 27TTO2 27TTOi 



qi = 



L., 



q-2 



l, 



2L :I 



or, in Cartesian coordinates, 
2nmi 



fci 



A-2 



2tt mi 



(A10) 



(All) 



V3aL x ' ZaLy/2 

The filling of the graphene Brillouin zone with discrete 
lattice momenta is illustrated on Fig. 1 in Subsection 
1.1. To obtain the complete system of eigenfunctions, it is 
sufficient to take m\ = . . . L x — 1, m 2 = . . . L y — 1. It is 
clear from (A10) that the Dirac points are only matched 
by discrete lattice momenta only if the lattice size L x is 
a multiple of 3 and L y is a multiple of 2. 

The eigenvalues and the eigenfunctions of the fermionic 
hopping matrices M a in the absence of interactions are 



'ipe ( s >£;<7) i 



A M (C, q,w) = l- e OTflr (1 - E c (q) At) , (A12) 

where "0C ( s j £; °) ar c the eigenfunctions of the one- 
particle Hamiltonian (Al). We note that for zero-energy 
states with Eq (q) = there is only one zero eigenvalue 
Am (C, q, to) of M a within the Brillouin zone w € [0, J^] , 
namely, at w = 0. Correspondingly, fermion propaga- 
tor in momentum space has only one pole. Thus the 
fermionic action (41) indeed describes a single fermion 
for each spin component a =J, \., and there is no fermion 
doubling problem. 
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Appendix B: Hopping expansion for the path 
integral representation of the partition function of 
the tight-binding model 

In this Appendix we consider the representation of 
the partition function (10) of the tight-binding model 
in terms of fermion worldlines in Euclidean space and 
discuss in more details the meaning of different approx- 
imations made in the derivation of the fermionic lattice 
action (41). With the help of the world- line representa- 
tion, we will prove that fluctuations of the electrostatic 
potential 0(s, £,r, z) cannot close the gap in the energy 
spectrum of the tight-binding model (9). In addition, wc 
will show that the lattice action (41) satisfies reflection 
positivity 42 , which is important for the self-consistency 
of the lattice regularization. 

Our starting point is the fermionic path integral with 
the weight which is a product of the factors (37) and the 

weights exp — fy CT (s,£) f] a (s,£) which come from 

V J 

the integral over the fermionic coherent states (33) in the 
decomposition of identity (32). The fermionic action can 
be then written as 

E E ( r ) M ° ( r > r ') *>* ( r ') . ( B1 ) 

CT=t,4, T,t' 



where 



M a (t, r ) = 5 TyT ' — 5 TyT i-\ exp (—h At ± 



(B2) 



and we have omitted the spatial coordinates s, £, treating 
the blocks of the fermion hopping matrix M„ (r, r') at 
fixed r, t' as operators which act on the space of single- 
particle wave functions if)(s,£). In order to account for 
the anti-periodic boundary conditions, the sign of the 
term proportional to <5 TjT '_i should be changed at, say, 
r = L T — 1. 

After integrating over the fermions, the partition func- 
tion (10) can be represented in the form (43). Consider 
now the logarithm of the determinant of M a in (43), 
assuming that it has the form (B2). Using the iden- 
tity logdet (M a ) = Tr log (M a ) and taking into account 
the special form of the fermion hopping matrix M a (r, r') 
with respect to r and r', we obtain 



logdet (M a ) = 



h T -\ 



logdet 1+ exp(-hAT±i(j)(T)) = 

t/At=0 



Y^^-Tr I II exp(-/iAr±^(r)) ] .(B3) 

it / At— 



The last expression can be interpreted as a sum over all 
possible configurations of a single fermionic world-line 
which wraps n times on the Euclidean torus with period 
(/cT)~ 1<)2 ' (ji , as illustrated on Fig. 15. For the world-line 



which originates from the site s, £ at time r and goes 
to the site s',£' at time r + Ar the weight is multiplied 
by the element [exp (— h At ± i(j> (t))] (s, £; s', £') of the 
single-particle transfer matrix. The factor (— 1) /n ac- 
counts for the Fermi statistics and compensates the over- 
counting of the world-line configurations due to ~ n dif- 
ferent ways to choose the starting time of the path. Obvi- 
ously, the full determinant det (M a ) = exp (logdet (M CT )) 
can be represented as a sum over any number N of 
fermionic world-lines, each coming with the weight (B3). 
As usual, an additional factor 1/Nl coming from the ex- 
pansion of the exponent then compensates for AH permu- 
tations of identical world-lines. 




FIG. 15: Fermionic worldlines which contribute to the par- 
tition function (10) with different approximations for the 
fermionic path integral. World-line 1 corresponds to the full 
single-particle transfer matrix exp ( — h At + i(f> (r)). Expan- 
sion up to the first order in Ar (38) allows only world-lines 
which hop once per Euclidean time interval At, such as world- 
line 2. If the integration over u in (38) is omitted, the world- 
lines can only hop at r = nAr, n £ Z, as the world-line 3. 

If wc expand the exponential exp (— h At ± i<f)) to 
the first order in Ar, as in (38), the worldlines are 
allowed to hop only once in the interval [t, t + At] . 
Each such hop changes the weight of the world-line by 
— h (s, £; s', £') At. For the tight-binding Hamiltonian 
(1), this means that only hops to nearest-neighbor sites 
on the hexagonal lattice are allowed with the weight 
k At. Time-like segments of worldlines contribute with 
the weight l±m At per time At depending on sublatticc 
index s. In addition, each world-line s (t) , £ (t) acquires 
the complex phase exp (-^i J drcj) (s (r) , £ (t) , t)) due to 
the presence of the electrostatic potential. 

Integration over u in (38) means that hops can happen 
at any time t + u At, u £ [0, 1] (see Fig. 15, world-line 2). 
Correspondingly, before the moment t+u At the fermion 
interacts with the electrostatic potential <j) (s, £,r, z = 0), 
and after that - with the potential (j> (s', r, z = 0). A 
simple estimate shows that if we neglect integration over 
u for two fermionic world-lines at distance r and replace 
the factors in (38) either by e iH^i.Tz=o) or e ^( a ',e',T,*=o) j 
this results in the correction of order O (nAre^Ar/r 2 ) 
to the path integral weight. Since we keep only the 
leading-order terms in kAt and e 2 a/r, we can discard 



25 



this term. However, if one uses some sort of improved 
action for the electromagnetic field which reproduces the 
Coulomb potential with higher accuracy, it might be nec- 
essary to keep the integration over u in (38). 

Let us now prove that the fluctuations of the elec- 
trostatic potential (f)(s,£,z — 0,r) cannot close the gap 
in the spectrum of the tight-binding Hamiltonian (9) at 
m^O. According to the discussion above, the logarithm 
of the determinant of the fermionic hopping matrix Mo- 
can be represented as the following sum over worldlines 
C: 



log dot (M a ) 



c 



(B4) 



where W [C] is the positive and real weight and e 4 *^ is 
the phase which includes both the integral of the electro- 
static potential along the world-line and the factor ( — 1)" 
in (B3). Due to the inequality | ^ aj| < l a »l one nas 



logdet(M CT )| <Y,W[C] 



(B5) 



By inverting the derivation above, it is easy to show that 
the latter sum corresponds to the bosonic partition func- 
tion 



^2W[C] =logdet(l -exp(-h/T)). 



(B6) 



which is finite if the single-particle Hamiltonian h has 
a gap. In Appendix A we have demonstrated that for 
to 7^ the spectrum of h indeed has a gap. From 
the inequality (B5) we see that log det (M ff ) remains fi- 
nite for any configuration of the electrostatic potential 
<j> (s, £, z, t), and thus the effective single-particle hamil- 
tonian h always has a gap. The finiteness of log det (M a ) 
at to ^ implies that det (M CT ) ^ and hence M a is 
invertible, which is crucial for the Hybrid Monte-Carlo 
algorithm. 

Finally, let us consider the reflection positivity of our 
lattice action. In our case, it is equivalent to the pos- 
itive definiteness of the single-particle transfer matrix 



exp (— hAr) . Its expansion 1 — hAr up to the first 
order in At is still positive-definite if the largest eigen- 
value of h does not exceed (At) -1 . From (A8) we see 
that this condition is equivalent to 



\Jrri 2 + 9k 2 At < 1. 



(B7) 



Appendix C: Current-current correlators in the free 
tight-binding model 



In this Appendix we consider the current-current cor- 
relators (59) and the corresponding spectral functions 
for the tight-binding model (9) without interactions with 
electromagnetic field. Our derivation is similar to that 
of 55 , but extends to the case of non-zero staggered po- 
tential to in (9). We start from the following general 
expression for the correlator of fermionic bilinear oper- 
ators J a — X) JA-.x,Y '$y m a theory with bilinear 
x.y 

Hamiltonian of the form H = Y] hx,Y wx ^ Y '- 

X,Y 



Z^Tr (,J A e- Tfl J B 



Tr [JA 



-rlh 



1 + e-P' 1 



Tr J B 



-08-T). 



1 + e-$ h 



Tr \j A 



-(f)-T)h 



1 



-ph 



JB 



1 + e-< ih 



(CI) 



where Z = Tr \ e~^ n J and (3 = T . The trace on the 

left-hand side is taken over the full Hilbert space of the 
theory, and on the right-hand side - over the one-particle 
Hilbert space (as in (35)). Applying the expression (CI) 
to the current-current correlator (59), we see that the 
contribution of the disconnected fermion diagrams (first 
summand on the right-hand side) is zero in this case. 
Evaluating the second term in the eigenbasis of the one- 
particle Hamiltonian (Al), we obtain 



e -(f3-T)H 



3V3L X L V l + e -Wr) i + e -/"Sc(«r 1 ' 

where the additional factor of two came from summation over spin indices and %,<&< (<Zi <?') is the- matrix element of 
the one-particle operator jf <a defined in (64) between the eigenstates (Al) of the one-particle Hamiltonian (35): 

ifc.ee (?) <?') = iK ( a > £; i) V'c W> £ + pb\ q') - $e (& ^ + pf>o) i>c' ( a > 6 ?') = 

= (wM a , 6 (q) Af PtC (q) e 1 ^ - mNp X (q) Af a , C (q) e"'^) L x L y 5 {q, q') . (C3) 
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Since our goal is to obtain the AC conductivity (57), we now contract the b and c indices of the correlator (C2) 
with the matrix Tb c introduced in (57): 

E T "c he c W, ?) he c (?> 1 1 ) = 3/2 Yl he c (?': s) jfc.c c (?> g')-i/2j2 he c (V* <?) ic,c c (?, <z') (C4) 

Explicit calculation yields the following expressions for the contracted matrix elements (C3) in (C4): 

J2jbXc(Q,9')=i/^Hq,q') (C-C) VE 2 (q) - m 2 (C5) 

5 



'EhC'cW.lihce M)= K 2&[q) ^ 2 (?) - 3CC'- 2 - (i? 2 (,) - m 2 ) CC Re (e™M£e-*«*jj (C6) 

Thus the squared matrix elements of the current operator are diagonal with respect to the lattice momenta and 
depend only on the product (,(,'. Denoting 



Yl Tbc he c (Y> 1) jcx c («r, ?') = I 



fas') C^C 

F++ (?)*(?,?') C = C 



with 



e 

b 

?2 



we can represent the current-current correlator (CI) in the following form: 



G (o) (t) = _JL_ V F+ - {q) C ° sh (2g (g) (T ~ P,2)) + F++ (g) 
{ 3 L x L y 2 cosh 2 (/3£;(g)/2) 



^ 2w V 4 / L ^ x 2/ , 4 cosh 2 (g) /2) 



(C7) 



F+ - (?) = 2V3^ ( g ) (9) + + ^ _ R ° ( e2i " (9) £ e " 2 " /Pb ) J _ E % " (C8) 



(C9) 



Comparing this expression with (60) and (61), we can obtain the AC conductivity en ' (w): 

^ H = E F + - (1) H»-U2 (?)) ^ tanh ( ^) + lllM. V (C10) 



One can see that the conductivity has a delta-function 
singularity at w = 0, which is a common feature of all 
ideal crystals due to the absence of scattering in an in- 
finite lattice without boundaries . In the limit m — > 
the non-singular part of the conductivity (first summand 
in (CIO)) reproduces the results of 55 . The finite part of 
a (w) is shown on Fig. 16. First, one can note a dis- 
tinct peak at w = 2 \J n 2 + m 2 , which corresponds to the 
singularity in the density of states associated with the 
saddle point at E = ^ n 2 + m 2 (point M in Fig. 1) in 
the dispersion relation. At m / and at sufficiently 
low temperatures, a second peak appears at w — 2m. 
In the intermediate frequency range, m, T -C w <C k, 
the conductivity approaches the value a = ao = 1 /4 
(l/4e 2 /ft = 7r/2e 2 //i in physical units). 



Let us now explicitly evaluate the expression (CIO) in 
the case when w is close to the threshold value w = 2m, 
so that only the momenta which are close to the Dirac 
points contribute to the summation over q in (C9) 
and (CIO). In this case we can approximate E (q) as 
yj m 2 + Vp fc 2 , where k is the momentum in Cartesian 
coordinates (A2). Assuming that the lattice size is suf- 
ficiently large, we can also replace summation over q 
in (C9), (CIO) (that is, summation over mi, m-i with 
q given by (A10)) by integration over k: L 27r L ^2 ~ 

* v q 

2 3v /^ a f d 2 k, where the factor of two accounts for two 

47T J 1 

Dirac points. A simple calculation shows that to the 
leading order in Sq = q — q^ the contributions of the 

terms proportional to Re ^e 2l ' a ^ e~ 2lqPb ^J have oppo- 
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T/k=0.02, m/K=0.1 
T/k=0.05, m/K=0.1 
T/k=0.1, m/K=0.1 
T/k=0.02, m/K=0 
T/k=0.05, m/K=0 
T/k=0.1, m/K=0 




DC conductivity of graphcnc oo = 1/4, which is cto = 
e 2 / (4h) = ne 2 / (2h) in SI units (see 10 for a discussion of 
the physical meaning of this result). From Fig. 16 one 
can see that when m = 0, this limiting value is reached at 
small frequencies w <^ k for temperatures T = 0.05 k and 
lower. On the other hand, the value of the conductivity 
at the threshold w = 2m is always two times larger than 
cr . 

Within the Dirac approximation one can also explicitly 
calculate the smeared low-frequency conductivity a 1 - " 1 
(introduced in (68)) in the absence of interactions: 



FIG. 16: AC conductivity on °> (w) (in units of ao — 
l/4e 2 /ft = 7r/2e 2 //i) for the non-interacting tight-binding 
model of graphene (9) with L x — L y — ► oo at different tem- 
peratures T and Dirac masses m. The inset shows the low- 
frequency region in a larger scale. 



site signs for q close to and and thus cancel in 
the sum over q in (C9) and (CIO) (see also' 5 ). Moreover, 
it is easy to show that the contribution of the last sum- 

mand in the expression (C8) for (q) is suppressed by 

an additional power of k a and thus can be also omitted 
in the Dirac approximation. With such simplifications, 
we can change the integration variable from k to E. In- 
tegrating out the delta-function in (CIO), we then arrive 
at the following expression for the finite part of the AC 
conductivity (CIO) for < (w — 2m) <^ k: 



-(o ) = ^ G (o) (/?/2) 



E 



F+_ (q) + F ++ (q) 



ttL x L v ^ 2 cosh z {/3E (q) /2) 



9 a 2 ac 2 /3 2 



4tt 2 



4 . 
73 lo S 



d 2 k- 



1 



rf 2 



cosh 2 (pE (k) j2) 

+ 00 



dEE 



cosh 2 (/3E/2) 

rn 

{ > ir 2 (1 + e m ) 



(C12) 



r(0) 



H«- i + 



4m" 



t anh ( — — 

4 



(Cll) 



If we take the limit j3 — > oo and m — > at finite w and 
then take w — 0, wc obtain the limiting value of the 



It is interesting to note that the value of in the limit 
m -> ct (0) = = 0.281 is quite close to the com- 

monly quoted value <jq — 1/4 and is also temperature- 
independent. 



* Electronic address: pavel.buividovich@physik.uni-regensbur; 
t Electronic address: polykarp@itep.ru 

1 A. K. Geim and K. S. Novoselov, Nature Materials 6, 183 
(2007), URL http://arxiv.org/abs/cond-mat/0702595. 

2 K. S. Novoselov, A. K. Geim, S. V. Morozov, 
D. Jiang, Y. Zhang, S. V. Dubonos, I. V. Grigorieva, 
and A. A. Firsov, Science 306, 666 (2004), URL 
http: //arxiv . org/abs/cond-mat/0410550. 

3 K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, 
M. I. Katsnelson, I. V. Grigorieva, S. V. Dubonos, 
and A. A. Firsov, Nature 438, 197 (2005), URL 
http : //arxiv . org/ abs/ cond-mat/0509330. 

4 Y. Zhang, Y. Tan, H. L. Stormer, and 
P. Kim, Nature 438, 201 (2005), URL 
http: //arxiv . org/ abs/ cond-mat/0509355. 

5 A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. 
Novoselov, and A. K. Geim, Rev. Mod. Phys. 81, 109 
(2009), URL http://arxiv.org/abs/0709.1163. 

6 G. W. Semenoff, Chiral symmetry breaking in 



.de graphene, Proceedings of the Nobel Symposium 
on Graphene and Quantum Matter (2011), URL 
http : //arxiv . org/abs/1 108 . 2945 . 

7 P. R. Wallace, Phys. Rev. 71, 622 (1947), URL 
http : //prola. aps . org/abstract/PR/v71/i9/p622_l . 

8 G. W. Semenoff, Phys. Rev. Lett. 53, 2449 (1984), URL 
http: //prl . aps . org/abstract/PRL/v53/i26/p2449_l . 

9 E. V. Gorbar, V. P. Gusynin, V. A. Miransky, and 
I. A. Shovkovy, Phys. Rev. B 66, 045108 (2002), URL 
http : //arxiv . org/abs/cond-mat/0202422 . 

10 M. I. Katsnelson, Eur. Phys. J. B 51, 157 (2006), URL 
http : //arxiv. org/abs/ cond-mat/0512337. 

11 J. E. Drut, T. A. Lahde, and E. Tolo, PoS Lattice2010, 
006 (2010), URL http://arxiv.org/abs/1011.0643. 

12 J. E. Drut and T. A. Lahde, Signatures of a 
gap in the conductivity of graphene (2010), URL 
http : //arxiv . org/abs/1005 . 5089 . 

13 G. Li, A. Luican, and E. Y. Andrei, Phys. Rev. Lett. 102, 
176804 (2009), URL http://arxiv.org/abs/0803.4016. 



28 



14 K. I. Bolotin, K. J. Sikes, J. Hone, H. L. Stormer, 
and P. Kim, Phys. Rev. Lett. 101, 096802 (2008), URL 
http : //arxiv . org/abs/0805 . 1830 . 

15 H. Leal and D. V. Khveshchenko, Nucl. Phys. B 687, 323 
(2004), URL http://arxiv.org/abs/cond-mat/0302164. 

16 O. V. Gamayun, E. V. Gorbar, and V. P. 



(2010), URL 



(2011), URL 



URL 



URL 



P. Blake, L. A. 
K. S. Novoselov, 
701 (2011), URL 



Gusynin, Phys. Rev. B 81, 075429 
http : //arxiv . org/abs/091 1 . 4878 . 

17 Y. Araki and T. Hatsuda, Phys. Rev. B 82, 121403 (2010), 
URL http://arxiv.org/abs/1003. 1769. 

18 Y. Araki, Ann.Phys. 326, 1408 
http : //arxiv . org/abs/1010 . 0847 . 

19 Y. Araki, Phys. Rev. B 85, 125436 (2012), 
http : //arxiv . org/abs/1201 . 1737 . 

20 D. T. Son, Phys. Rev. B 75, 235423 (2007), 
http : //arxiv . org/abs/cond-mat/0701501 . 

21 D. C. Elias, R. V. Gorbachev, A. S. Mayorov, 
S. V. Morozov, A. A. Zhukov, 
Ponomarenko, I. V. Grigorieva, 
F. Guinea, et al., Nature Phys. 7, 
http : //arxiv . org/abs/1 104 . 1396 . 

22 A. S. Mayorov, D. C. Elias, I. S. Mukhin, S. V. Moro- 
zov, L. A. Ponomarenko, K. S. Novoselov, A. K. Geim, 
and R. V. Gorbachev, How close can one approach the 
Dirac point in graphene experimentally'? (2012), URL 
http : //arxiv . org/abs/1206 . 3848 . 

23 R. Shankar, Rev. Mod. Phys. 66, 129192 (1994), URL 
http : //rmp . aps . org/abstract/RMP/v66/il/pl29_l . 

24 S. A. Jafari, Eur.Phys.J.B 68, 537 (2009), URL 
http : //arxiv . org/abs/0809 . 1109 . 

25 J. E. Drut and D. T. Son, Phys. Rev. B 77, 075115 (2008), 
URL http://arxiv.org/abs/0710. 1315. 

26 J. E. Drut and T. A. Lahde, Phys. Rev. Lett. 102, 026802 
(2009), URL http://arxiv.org/abs/0807.0834. 

27 J. E. Drut and T. A. Lahde, Phys. Rev. B 79, 
(2009), URL http://arxiv.org/abs/0901.0584. 

28 J. E. Drut and T. A. Lahde, Phys. Rev. B 79, 
(2009), URL http://arxiv.org/abs/0905.1320. 

29 J. E. Drut, T. A. Lahde, and L. Suoranta, First- 
order chiral transition in the compact lattice theory of 
graphene and the case for improved actions (2010), URL 
http : //arxiv . org/abs/1002 . 1273 . 

30 J. E. Drut and T. A. Lahde, PoS Lattice2011, 074 (2011), 
URL http : //arxiv . org/abs/1 1 1 1 . 0929 . 

31 S. Hands and C. Strouthos, Phys. Rev. B 78, 165423 
(2008), URL http://arxiv.org/abs/0806.4877. 

32 W. Armour, S. Hands, and C. Strouthos, Phys. Rev. B 81, 
125105 (2010), URL http://arxiv.org/abs/0910.5646. 

33 W. Armour, S. Hands, and C. Strouthos, Phys. Rev. B 84, 
075123 (2011), URL http://arxiv.org/abs/1105.1043. 

34 P. V. Buividovich, E. V. Luschevskaya, O. V. Pavlovsky, 
M. I. Polikarpov, and M. V. Ulybyshev, Numeri- 
cal study of the conductivity of graphene monolayer 
within the effective field theory approach (2012), URL 
http : //arxiv . org/abs/1204 . 0921 . 

35 J. Giedt, A. Skinner, and S. Nayak, Phys. Rev. B 83, 
045420 (2011), URL http://arxiv.org/abs/0911.4316. 

36 G. Metalidis, Ph.D. thesis, University of 
Halle- Wittenberg, Germany (2007), URL 
http : //sundoc . bibliothek .uni-halle . de/diss-online/07/07H©i£@iton 

37 R. Brower, C. Rebbi, and D. Schaich, PoS LAT2011, 056 
(2012), URL http://arxiv.org/abs/1204.5424. 

38 R. C. Brower, C. Rebbi, and D. Schaich, Hybrid Monte 081409 (2011), URL http: //arxiv . org/abs/1006 . 2426. 
Carlo simulation of graphene on the hexagonal lattice 62 B. B. Beard and U. J. Wiese, Phys. Rev. Lett. 77, 5130 



165425 



241405 



(2011), URL http://arxiv.org/abs/1101.5131. 

39 K. Nomura, S. Ryu, and D. H. Lee, Phys. Rev. Lett. 103, 
216801 (2009), URL http://arxiv.org/abs/0906.0159. 

40 C. Burden and A. N. Burkitt, Eur. 
Phys. Lett. 3, 545 (1987), URL 
http : //dx . doi . org/10 . 1209/0295-5075/3/5/006 . 

41 M. Creutz, Phys. Lett. B 649, 230 (2007), URL 
http://arxiv.org/abs/hep-lat/0701018. 

42 I. Montvay and G. Muenster, Quantum fields on a lattice 
(Cambridge University Press, 1994). 

43 K. V. Zakharchenko, A. Fasolino, J. H. Los, and M. I. 
Katsnelson, J.Phys.Cond.Mat. 23, 202202 (2011), URL 
http://arxiv.org/abs/1104. 1130. 

44 R. E. Peierls, Z. Phys. 80, 763 (1933), URL 
http : //dx . doi . org/10 . 1007/BF01342591 . 

45 D. R. Hofstadter, Phys. Rev. B 14, 2239 (1976), URL 
http : //prb . aps . org/abstract/PRB/vl4/i6/p2239_l . 

46 D. Chakrabarti, S. Hands, and A. Rago, JHEP 0906, 060 
(2009), URL http://arxiv.org/abs/0904.1310. 

47 T. O. Wehling, E. §a§ioglu, C. Friedrich, A. I. Lichtenstein, 
M. I. Katsnelson, and S. Bliigel, Phys. Rev. Lett. 106, 
236805 (2011), URL http://arxiv.org/abs/1101.4007. 

48 T. DeGrand and C. DeTar, Lattice methods for quantum 
chromodynamics (World Scientific, 2006). 

49 S. Gottlieb, W. Liu, D. Toussaint, R. L. Renken, and 
R. L. Sugar, Phys. Rev. D 35, 2531 (1987), URL 
http : //prd. aps . org/abstract/PRD/v35/i8/p2531_l . 

50 J. C. Sexton and D. H. Weingarten, 
Nucl. Phys. B 380, 665 (1992), URL 
http : //dx . doi . org/10 . 1016/0550-3213 (92) 90263-B. 

51 L. Del Debbio and S. Hands, Phys. Lett. B 373, 171 
(1996), URL http://arxiv.org/abs/hep-lat/9512013. 

52 M. Liischer, Comput.Phys.Commun. 79, 100 (1994), URL 
http : //arxiv . org/abs/hep-lat /9309020 . 

53 C. Bernard, S. Hashimoto, D. B. Leinweber, P. Lep- 
age, E. Pallante, S. R. Sharpe, and H. Wittig, Nucl. 
Phys. B - Proc. Suppl. 119, 170 (2003), URL 
http : //arxiv . org/abs/hep-lat/0209086 . 

54 C. Alexandrou, V. Drach, K. Hadjiyiannakou, 
K. Jansen, G. Koutsou, A. Strelchenko, and 
A. Vaquero, PoS Lattice2012, 184 (2012), URL 
http : //arxiv . org/abs/1211 . 0126 . 

55 T. Stauber, N. M. R. Peres, and A. K. 
Geim, Phys. Rev. B 78, 085432 (2008), URL 
http : //arxiv . org/abs/0803 . 1802 . 

56 L. P. Kadanoff and P. C. Martin, 
Ann. Phys. 24, 419 (1963), URL 
http://dx.doi.org/10. 1016/0003-4916(63)90078-2. 

57 A. Hosoya, M. Sakagami, and M. Takao, 
Ann.Phys. 154, 229 (1984), URL 
http : //dx . doi . org/10 . 1016/0003-4916 (84) 90144- 1 . 

58 M. Asakawa, T. Hatsuda, and Y. Nakahara, 
Prog. Part. Nucl. Phys. 46, 459 (2001), URL 
http://arxiv.org/abs/hep-lat/0011040. 

59 G. Aarts, C. Allton, J. Foley, S. Hands, and 
S. Kim, Phys. Rev. Lett. 99, 022002 (2007), URL 
http : //arxiv . org/abs/hep- lat /0703008 . 

60 B. Rosenstein, M. Lewkowicz, and T. Maniv, 
Chiral anomaly and strength of the electron- 
interaction in graphene (2012), URL 



http : //arxiv . org/abs/ 1210 . 3345 . 
61 M. van Schilfgaarde and M. I. Katsnelson, Phys. Rev. B 83, 



(1996), URL http : //arxiv . org/abs/cond-mat/9602164. http : //dx . doi . org/10 . 1016/0003-4916 (92) 90288-W. 

E. Farhi and S. Gutmann, Ann. 64 We thank Prof. Dr. M. I. Katsnelson for pointing this 

Phys. 213, 182 (1992), URL us. 



